\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}
\usepackage[
  colorlinks=true,
  urlcolor=BlueViolet,
  citecolor=BlueViolet,
  filecolor=BlueViolet,
  linkcolor=BlueViolet
]{hyperref}


\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 nearly 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 in random landscapes on
  their marginality, and apply it in three isotropic settings with
  qualitatively different structure: in the spherical spin-glasses, where the
  energy is Gaussian and its Hessian is GOE; in a multispherical spin glasses,
  which are Gaussian but non-GOE; and in spherical 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 statistical inference problems. The
behavior of these systems is best understood when equilibrium or optimal
solutions are studied and weighted 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 mean-field settings, it was long thought that physical and many algorithmic
dynamics would get stuck at a specific energy level, called the threshold
energy. The threshold energy is the energy level at which level sets of the
landscape transition from containing mostly saddle points to containing mostly
minima. At this threshold, the level set contains mostly \emph{marginal
minima}, or minima that have a pseudogap in the spectrum of their Hessian.

However, recent work found that the threshold energy is not important even for
simple gradient descent dynamics \cite{Folena_2020_Rethinking, Folena_2023_On}.
Depending on the initial condition of the system and the nature of the
dynamics, the energy reached can be above or below the threshold energy, while
in some models the threshold energy is completely inaccessible to any dynamics
\cite{Kent-Dobias_2023_How}. Though it is still not known how to predict the
energy level that many simple algorithms will reach, the results all share one
commonality: the minima found are still marginal, despite being in the minority
compared to stiff minima or saddle points. This ubiquity of behavior suggests
that the distribution of marginal minima can bound out-of-equilibrium dynamical
behavior.

It is not straightforward to condition on the marginality of minima using the
traditional methods for analyzing the distribution of minima in rugged
landscapes. Using the method of a Legendre transformation of the Parisi
parameter corresponding to a set of real replicas, one can force the result to
correspond with marginal minima by tuning the value of that parameter, but this
results in only the threshold energy and cannot characterize marginal minima at
energies where they are a minority \cite{Monasson_1995_Structural}.

The alternative approach, used to great success in the spherical spin glasses, is to
start by making a detailed 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 spin glasses because it is straightforward to implement.
First, the shape of the Hessian's spectrum is independent of energy and even
whether one sits at a stationary point or not. This is a property of models
whose energy is a Gaussian random variable \cite{Bray_2007_Statistics}.
Furthermore, a natural parameter in the analysis of these models linearly
shifts the spectrum of the Hessian. Therefore, tuning this parameter to a
specific constant value allows one to require that the Hessian spectrum have a
pseudogap, and therefore that the associated minima be marginal. Unfortunately
this strategy is less straightforward to generalize to other models. Many models of interest,
especially in inference problems, have Hessian statistics that are poorly
understood. This is especially true for the statistics of the Hessian
conditioned to lie at stationary points, which is necessary to understand in
models whose energy is non-Gaussian.

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 integration measure to select only configurations that result
in a certain value of the smallest eigenvalue of a matrix. By requiring that
the smallest eigenvalue of the Hessian at stationary points be zero, and
further looking for a sign that the zero eigenvalue lies at the edge of a
continuous spectrum, we enforce the condition that the spectrum has a
pseudogap, and is therefore marginal. We demonstrate the method on the
spherical spin glasses, where it is unnecessary but instructive, and on extensions of
the spherical models with non-GOE Hessians where the technique is more useful.
In a companion paper, we compare the marginal complexity with the performance
of gradient descent and approximate message passing \cite{Kent-Dobias_2024_Algorithm-independent}.

In Section \ref{sec:eigenvalue} we introduce the technique for conditioning on
the smallest eigenvalue and how to extend it to further condition on the
presence of a pseudogap. We provide a simple but illustrative example using a
shifted GOE matrix. In Section \ref{sec:marginal.complexity} we apply this
technique to the problem of characterizing marginal minima in random
landscapes. The following Section \ref{sec:examples} gives several examples of
the marginal complexity applied to specific models of increasing difficulty.
Finally, Section \ref{sec:conclusion} summarizes this work and suggests
necessary extensions.

\section{Conditioning on the smallest eigenvalue}
\label{sec:eigenvalue}

In this section, we introduce a general method for conditioning a measure on
the smallest eigenvalue of some matrix that depends on it. In Section
\ref{sec:shifted.GOE} we show how this works in perhaps the simplest example of
GOE random matrices with a shifted diagonal. In the final subsection we
describe how to extend this method to condition on the presence of a pseudogap
at the bottom on the spectrum.

\subsection{The general method}

An arbitrary function $g$ of the minimum eigenvalue of a matrix $A$ can be
expressed as
\begin{equation} \label{eq:λmin}
  g(\lambda_\textrm{min}(A))
  =\lim_{\beta\to\infty}\int
    \frac{d\mathbf s\,\delta(N-\|\mathbf s\|^2)e^{-\beta\mathbf s^TA\mathbf s}}
      {\int d\mathbf s'\,\delta(N-\|\mathbf s'\|^2)e^{-\beta\mathbf s'^TA\mathbf s'}}
    g\left(\frac{\mathbf s^TA\mathbf s}N\right)
\end{equation}
In the limit of large $\beta$, each integral concentrates among vectors
$\mathbf s$ in the eigenspace of $A$ corresponding to the smallest eigenvalue
of $A$. This produces
\begin{equation}
  \begin{aligned}
    &\lim_{\beta\to\infty}\int\frac{
      d\mathbf s\,\delta(N-\|\mathbf s\|^2)e^{-\beta\mathbf s^TA\mathbf s}
    }{
      \int d\mathbf s'\,\delta(N-\|\mathbf s'\|^2)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\|^2)\mathbb 1_{\operatorname{ker}(A-\lambda_\mathrm{min}(A)I)}(\mathbf s)
    }{
      \int d\mathbf s'\,\delta(N-\|\mathbf s'\|^2)\mathbb 1_{\operatorname{ker}(A-\lambda_\mathrm{min}(A)I)}(\mathbf s')}g\left(\frac{\mathbf s^TA\mathbf s}N\right) \\
    &=g(\lambda_\mathrm{min}(A))
    \frac{\int d\mathbf s\,\delta(N-\|\mathbf s\|^2)\mathbb 1_{\operatorname{ker}(A-\lambda_\mathrm{min}(A)I)}(\mathbf s)}{\int d\mathbf s'\,\delta(N-\|\mathbf s'\|^2)\mathbb 1_{\operatorname{ker}(A-\lambda_\mathrm{min}(A)I)}(\mathbf s')} \\
    &=g(\lambda_\mathrm{min}(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} in the context of random landscapes. 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}
\label{sec: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{min}(B+\mu I)=\lambda^*}
  =\overline{\delta\big(N\lambda^*-N\lambda_\mathrm{min}(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{min}$ defined in \eqref{eq:λmin}, we have
\begin{widetext}
\begin{equation}
  e^{NG_{\lambda^*}(\mu)}
  =\overline{
    \lim_{\beta\to\infty}\int\frac{d\mathbf s\,\delta(N-\|\mathbf s\|^2)e^{-\beta\mathbf s^T(B+\mu I)\mathbf s}}
    {\int d\mathbf s'\,\delta(N-\|\mathbf s'\|^2)e^{-\beta\mathbf s'^T(B+\mu I)\mathbf s'}}\,\delta\big(N\lambda^*-\mathbf s^T(B+\mu I)\mathbf s\big)
  }
\end{equation}
Using replicas to treat the denominator ($x^{-1}=\lim_{m\to0}x^{m-1}$)
and transforming the $\delta$-function to its Fourier
representation, we have
\begin{equation}
  e^{NG_{\lambda^*}(\mu)}
  =\overline{\lim_{\beta\to\infty}\lim_{m\to0}\int d\hat\lambda\prod_{\alpha=1}^m\left[d\mathbf s^\alpha\,\delta(N-\|\mathbf s^\alpha\|^2)\right]
  \exp\left\{-\beta\sum_{\alpha=1}^m(\mathbf s^\alpha)^T(B+\mu I)\mathbf s^\alpha+\hat\lambda\left[N\lambda^*-(\mathbf s^1)^T(B+\mu I)\mathbf s^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_{m\to0}\int d\hat\lambda\prod_{\alpha=1}^m\left[d\mathbf s^\alpha\,\delta(N-\|\mathbf s^\alpha\|^2)\right] \\
  &\hspace{10em}\exp\left\{N\left[\hat\lambda(\lambda^*-\mu)-m\beta\mu\right]+\frac{\sigma^2}{N}\left[\beta^2\sum_{\alpha\gamma}^m(\mathbf s^\alpha\cdot\mathbf s^\gamma)^2
    +2\beta\hat\lambda\sum_\alpha^m(\mathbf s^\alpha\cdot\mathbf s^1)^2
    +\hat\lambda^2N^2
  \right]\right\}
  \end{aligned}
\end{equation}
\end{widetext}
We make the Hubbard--Stratonovich transformation to the matrix field
$Q_{ab}=\frac1N\mathbf s_a^T\mathbf s_b$. This gives
\begin{equation}
  e^{NG_{\lambda^*}(\mu)}
  =\lim_{\beta\to\infty}\lim_{m\to0}\int d\hat\lambda\,dQ\,
  e^{N\mathcal U_\mathrm{GOE}(\hat\lambda,Q\mid\mu,\lambda^*,\beta)}
\end{equation}
where the effective action is given by
\begin{equation}
  \begin{aligned}
    &\mathcal U_\textrm{GOE}(\hat\lambda, Q\mid\lambda^*,\mu,\beta)
    =\hat\lambda(\lambda^*-\mu)-m\beta\mu \\
    &+\sigma^2\left[\beta^2\sum_{ab}^mQ_{ab}^2
        +2\beta\hat\lambda\sum_a^mQ_{1a}^2
      +\hat\lambda^2
    \right]+\frac12\log\det Q
  \end{aligned}
\end{equation}
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} \label{eq:Q.structure}
  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&1
  \end{bmatrix}
\end{equation}
with $\sum_{ab}Q_{ab}^2=m+2(m-1)\tilde q_0^2+(m-1)(m-2)q_0^2$, $\sum_aQ_{1a}^2=1+(m-1)\tilde q_0^2$,
and
\begin{equation}
  \log\det Q=(m-2)\log(1-q_0)+\log(1+(m-2)q_0-(m-1)\tilde q_0^2)
\end{equation}
Inserting these expressions and taking the limit of $m$ 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 U_\textrm{GOE}(\hat\lambda,q_0,\tilde q_0\mid\mu,\lambda^*,\beta)}
\end{equation}
with the effective action
\begin{equation}
  \begin{aligned}
    &\mathcal U_\mathrm{GOE}(\hat\lambda,q_0,\tilde q_0\mid\mu,\lambda^*,\beta) \\
    &\quad=\hat\lambda(\lambda^*-\mu)+\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}
  \label{eq:q0.limit}
  q_0&=1-y\beta^{-1}-z\beta^{-2}+O(\beta^{-3})
  \\
  \label{eq:q0t.limit}
  \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 U_\textrm{GOE}(\hat\lambda,y,\Delta z\mid\mu,\lambda^*,\infty)
    &=\hat\lambda(\lambda^*-\mu)
    +\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$, 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} \label{eq:goe.large.dev}
  \begin{aligned}
    &G_{\lambda^*}(\mu)
    =\mathop{\textrm{extremum}}_{y,\Delta z,\hat\lambda}
    \mathcal U_\mathrm{GOE}(y,\Delta z,\hat\lambda\mid\mu,\lambda^*,\infty) \\
    &=-\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 $\lambda^*=0$. For $\mu<2\sigma$ $G_{0}(\mu)$ has an
imaginary part. 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}
  \hspace{1.3em}
  \includegraphics{figs/spectrum_less.pdf}
  \hspace{-2em}
  \includegraphics{figs/spectrum_eq.pdf}
  \hspace{-2em}
  \includegraphics{figs/spectrum_more.pdf}
  \\
  \includegraphics{figs/large_deviation.pdf}
  \caption{
    The large deviation function $G_0(\mu)$ defined in
    \eqref{eq:large.dev} as a function of the shift $\mu$ to the
    GOE diagonal. $G_0(2\sigma)=0$, while for
    $\mu>2\sigma$ it is negative and for $\mu<2\sigma$ it gains an
    imaginary part. The top panels show schematically what happens to the
    spectral density in each of these regimes. For $\mu<2\sigma$, an $N^2$
    large deviation would be required to fix the smallest eigenvalue to zero
    and the calculation breaks down, leading to the imaginary part. For
    $\mu>2\sigma$ the spectrum can satisfy the constraint on the smallest
    eigenvalue by isolating a single eigenvalue at zero at the cost of an
    order-$N$ large deviation. At the transition point $\mu=2\sigma$ the
    spectrum is pseudogapped.
  } \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
touching 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 example problem of section \ref{sec:shifted.GOE}, 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}
\label{sec:marginal.complexity}

The methods of the previous section can be used in diverse settings. However,
we are interested in applying them to study stationary points in random
landscapes whose Hessian spectrum has a pseudogap -- that is, that are
marginal. In Section \ref{sec:marginal.kac-rice} we define the marginal
complexity using the tools of the previous section. In Section
\ref{sec:general.features} we review several general features in a physicists'
approach to computing the marginal complexity. In Section
\ref{sec:superspace_kac-rice} we introduce a representation of the marginal
complexity in terms of an integral over a superspace, which condenses the
notation and the resulting calculation and which we will use in one of our
examples in the next section.

\subsection{Marginal complexity from Kac--Rice}
\label{sec:marginal.kac-rice}

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 x)=0$ for $\mathbf x\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 the Lagrangian
\begin{equation}
  L(\mathbf x,\pmb\omega)=H(\mathbf x)+\sum_{i=1}^r\omega_ig_i(\mathbf x)
\end{equation}
with respect to $\mathbf x$ and $\pmb\omega=\{\omega_1,\ldots,\omega_r\}$. The corresponding gradient and Hessian for the problem are
\begin{align}
  \nabla H(\mathbf x,\pmb\omega)
  &=\partial L(\mathbf x,\pmb\omega)
  =\partial H(\mathbf x)+\sum_{i=1}^r\omega_i\partial g_i(\mathbf x)
  \\
  \operatorname{Hess}H(\mathbf x,\pmb\omega)
  &=\partial\partial L(\mathbf x,\pmb\omega)
  =\partial\partial H(\mathbf x)+\sum_{i=1}^r\omega_i\partial\partial g_i(\mathbf x)
\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} \label{eq:kac-rice.measure}
  d\nu_H(\mathbf x,\pmb\omega)=d\mathbf x\,d\pmb\omega\,\delta\big(\nabla H(\mathbf x,\pmb\omega)\big)\,\delta\big(\mathbf g(\mathbf x)\big)\,\big|\det\operatorname{Hess}H(\mathbf x,\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 as
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 and spectrum trace,
\begin{equation} \label{eq:kac-rice.measure.2}
  \begin{aligned}
    &d\nu_H(\mathbf x,\pmb\omega\mid E,\mu) \\
    &\quad=d\nu_H(\mathbf x,\pmb\omega)\,
    \delta\big(NE-H(\mathbf x)\big)
    \,\delta\big(N\mu-\operatorname{Tr}\operatorname{Hess}H(\mathbf x,\pmb\omega)\big)
  \end{aligned}
\end{equation}
We further want to control the value of the minimum eigenvalue of the Hessian at the stationary points. Using the method introduced above, we can write the number of stationary points with energy $E$, Hessian trace $\mu$, and smallest eigenvalue $\lambda^*$ as
\begin{widetext}
\begin{equation}
  \begin{aligned}
    &\mathcal N_H(E,\mu,\lambda^*)
    =\int d\nu_H(\mathbf x,\pmb\omega\mid E,\mu)\,\delta\big(N\lambda^*-\lambda_\mathrm{min}(\operatorname{Hess}H(\mathbf x,\pmb\omega))\big) \\
    &=\lim_{\beta\to\infty}\int d\nu_H(\mathbf x,\pmb\omega\mid E,\mu)
    \frac{d\mathbf s\,\delta(N-\|\mathbf s\|^2)\delta(\mathbf s^T\partial\mathbf g(\mathbf x))e^{-\beta\mathbf s^T\operatorname{Hess}H(\mathbf x,\pmb\omega)\mathbf s}}
    {\int d\mathbf s'\,\delta(N-\|\mathbf s'\|^2)\delta(\mathbf s'^T\partial\mathbf g(\mathbf x))e^{-\beta\mathbf s'^T\operatorname{Hess}H(\mathbf x,\pmb\omega)\mathbf s'}}
    \delta\big(N\lambda^*-\mathbf s^T\operatorname{Hess}H(\mathbf x,\pmb\omega)\mathbf s\big)
  \end{aligned}
\end{equation}
where the $\delta$-functions
\begin{equation}
  \delta(\mathbf s^T\partial\mathbf g(\mathbf x))
  =\prod_{s=1}^r\delta(\mathbf s^T\partial g_i(\mathbf x))
\end{equation}
ensure that the integrals are constrained to the tangent space of the
configuration manifold at the point $\mathbf x$. This likewise allows us to
define the complexity of points with a specific energy, stability, and minimum 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 ($x^{-1}=\lim_{m\to-1}x^m$). This leads to the expression
\begin{equation} \label{eq:min.complexity.expanded}
  \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\nu_H(\mathbf x_a,\pmb\omega_a\mid E,\mu)\,\delta\big(N\lambda^*-(\mathbf s_a^1)^T\operatorname{Hess}H(\mathbf x_a,\pmb\omega_a)\mathbf s_a^1\big)\\
    &\hspace{12em}\times\lim_{m_a\to0}
      \left(\prod_{\alpha=1}^{m_a} d\mathbf s_a^\alpha
      \,\delta\big(N-\|\mathbf s_a^\alpha\|^2\big)
      \,\delta\big((\mathbf s_a^\alpha)^T\partial\mathbf g(\mathbf x_a)\big)
      \,e^{-\beta(\mathbf s_a^\alpha)^T\operatorname{Hess}H(\mathbf x_a,\pmb\omega_a)\mathbf s_a^\alpha}\right)
    \Bigg]
  \end{aligned}
\end{equation}
\end{widetext}
for the complexity of stationary points of a given energy, trace, and smallest eigenvalue.

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 minimum 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=\mu_\text{m}(E)$,
\begin{equation}
  \Sigma_\text{m}(E)
  =\Sigma_0(E,\mu_\text m(E))
\end{equation}

\subsection{General features of saddle point computation}
\label{sec:general.features}

Several elements of the computation of the marginal complexity, and indeed the
ordinary dominant complexity, follow from the formulae of the above section in
the same way. The physicists' approach to this problem seeks to convert all of
the Kac--Rice measure defined in \eqref{eq:kac-rice.measure} and
\eqref{eq:kac-rice.measure.2} into elements of an exponential integral over
configuration space. To begin with, all Dirac $\delta$ functions are
expressed using their Fourier representation, with
\begin{align}
  \label{eq:delta.grad}
  &\delta\big(\nabla H(\mathbf x_a,\pmb\omega_a)\big)
    =\int\frac{d\hat{\mathbf x}_a}{(2\pi)^N}e^{i\hat{\mathbf x}_a^T\nabla H(\mathbf x_a,\pmb\omega_a)} \\
    \label{eq:delta.energy}
  &\delta\big(NE-H(\mathbf x_a)\big)
    =\int\frac{d\hat\beta_a}{2\pi}e^{\hat\beta_a(NE-H(\mathbf x_a))} \\
  \label{eq:delta.eigen}
  &\begin{aligned}
    &\delta\big(N\lambda^*-(\mathbf s_a^1)^T\operatorname{Hess}H(\mathbf x_a,\pmb\omega)\mathbf s_a^1\big) \\
    &\qquad\qquad\qquad=\int\frac{d\hat\lambda_a}{2\pi}e^{\hat\lambda_a(N\lambda^*-(\mathbf s_a^1)^T\operatorname{Hess}H(\mathbf x_a,\pmb\omega)\mathbf s_a^1)}
  \end{aligned}
\end{align}
To do this we have introduced auxiliary fields $\hat{\mathbf x}_a$,
$\hat\beta_a$, and $\hat\lambda_a$. Since the permutation symmetry of vector
elements is preserved in \textsc{rsb} order, the order parameters $\hat\beta$
and $\hat\lambda$ will quickly lose their indices, since they will ubiquitously
be constant over the replicas at the eventual saddle point solution.

We would like to make a similar treatment of the determinant of the Hessian
that appears in \eqref{eq:kac-rice.measure}. The standard approach is to drop
the absolute value function around the determinant. This can potentially lead
to severe problems with the complexity. However, it is a justified step when
the parameters of the problem, i.e., $E$, $\mu$, and $\lambda^*$ put us in a
regime where the exponential majority of stationary points have the same index.
This is true for maxima and minima, and for saddle points whose spectra have a
strictly positive bulk with a fixed number of negative outliers. It is in
particular a safe operation for this problem of marginal minima, which lie
right at the edge of disaster. Dropping the absolute value sign allows us to
write
\begin{equation} \label{eq:determinant}
  \det\operatorname{Hess}H(\mathbf x_a, \pmb\omega_a)
  =\int d\bar{\pmb\eta}_a\,d\pmb\eta_a\,e^{-\bar{\pmb\eta}_a^T\operatorname{Hess}H(\mathbf x_a,\pmb\omega)\pmb\eta_a}
\end{equation}
for $N$-dimensional Grassmann variables $\bar{\pmb\eta}_a$ and $\pmb\eta_a$. For
the spherical models this step is unnecessary, since there are other ways to
treat the determinant keeping the absolute value signs, as in previous works
\cite{Folena_2020_Rethinking, Kent-Dobias_2023_How}. However, other of
our examples are for models where the same techniques are impossible.

For the cases studied here, fixing the trace results in a relationship
between $\mu$ and the Lagrange multipliers enforcing the constraints. This is
because the trace of $\partial\partial H$ is typically an order of $N$ smaller
than that of the constraint functions $\partial\partial g_i$. The result is that
\begin{equation}
  \mu
  =\frac1N\operatorname{Tr}\operatorname{Hess}H(\mathbf x)
  =\frac1N\sum_{i=1}^r\omega_i\partial\partial g_i(\mathbf x)
  +O(N^{-1})
\end{equation}
In particular, here we study only cases with quadratic $g_i$, which results in
an expression relating $\mu$ and the $\omega_i$ that is independent of $\mathbf
x$. Since $H$ contains the disorder of the problem, this simplification means
that the effect of fixing the trace is independent of the disorder and only
depends on properties of the constraint manifold.

\subsection{Superspace representation}
\label{sec:superspace_kac-rice}

The ordinary Kac--Rice calculation involves many moving parts, and this method
for incorporating marginality adds even more. It is therefore convenient to
introduce compact and simplifying notation through a superspace representation.
The use of superspace in the Kac--Rice calculation is well established, as well
as the deep connections with BRST symmetry that is implied.
Appendix~\ref{sec:superspace} introduces the notation and methods of
superspace. Here we describe how it can be used to simplify the complexity
calculation in the marginal case.

We consider the $\mathbb R^{N|4}$ superspace whose Grassmann indices are
$\bar\theta_1,\theta_1,\bar\theta_2,\theta_2$. Consider the supervector defined
by
\begin{equation}
  \pmb\phi_a^\alpha(1,2)
  =\mathbf x_a
  +\bar\theta_1\pmb\eta_a+\bar{\pmb\eta}_a\theta_1
  +i\hat{\mathbf x}_a\bar\theta_1\theta_1
  +\mathbf s_a^\alpha(\bar\theta_1\theta_2+\bar\theta_2\theta_1)
\end{equation}
Note that this supervector does not span the whole superspace: only a couple
terms from the $\bar\theta_2,\theta_2$ sector are present, since the rest are
unnecessary for our representation. With this supervector so defined, the
replicated count of stationary points with energy $E$, trace $\mu$, and
smallest eigenvalue $\lambda^*$ can be written as
\begin{widetext}
\begin{equation}
  \begin{aligned}
    \mathcal N_H(E,\mu,\lambda^*)^n
    &=\lim_{\beta\to\infty}\int\prod_{a=1}^nd\pmb\omega_a\lim_{m_a\to0}\prod_{\alpha=1}^{m_a}d\pmb\phi_a^\alpha
    \exp\left\{
      \delta^{\alpha1}N(\hat\beta_aE+\hat\lambda_a\lambda^*)
      +\int d1\,d2\,B_a^\alpha(1,2)L(\pmb\phi_a^\alpha(1,2),\pmb\omega_a)
    \right\}
  \end{aligned}
\end{equation}
Here we have also defined the operator
\begin{equation}
  B_a^\alpha(1,2)=\delta^{\alpha1}\bar\theta_2\theta_2
        (1-\hat\beta_a\bar\theta_1\theta_1)
        -\delta^{\alpha1}\hat\lambda_a-\beta
\end{equation}
which encodes various aspects of the complexity problem, and the measures
\begin{align}
  d\pmb\phi_a^\alpha
  &=\left[
    d\mathbf x_a\,\delta\big(\mathbf g(\mathbf x_a)\big)\,
    \frac{d\hat{\mathbf x}_a}{(2\pi)^N}\,
    d\pmb\eta_a\,d\bar{\pmb\eta}_a\,
    \delta^{\alpha1}+(1-\delta^{\alpha1})
  \right]\,
  d\mathbf s_a^\alpha\,\delta(\|\mathbf s_a^\alpha\|^2-N)\,
  \delta\big((\mathbf s_a^\alpha)^T\partial\mathbf g(\mathbf x_a)\big)
  \\
  d\pmb\omega_a&=\prod_{i=1}^rd\omega_{ai}\,\delta\big(N\mu-\omega_{ai}\partial\partial g_i(\mathbf x_a)\big)
\end{align}
that collect the individual measures of the various fields embedded in the superfield.
\end{widetext}
With this way of writing the replicated count, the problem of marginal
complexity temporarily takes the schematic form of an equilibrium calculation
with configurations $\pmb\phi$, inverse temperature $B$, and energy $L$. This
makes the intermediate pieces of the calculation dramatically simpler. Of
course the complexity of the underlying problem is not banished: near the end
of the calculation, terms involving the superspace must be expanded.

\section{Examples}
\label{sec:examples}

In this section we present analysis of marginal complexity in three random
landscapes. In Section \ref{sec:ex.spherical} we apply the methods described
above to the spherical spin glasses, which reveals some general aspects of the
calculation. Since the spherical spin glasses are Gaussian and have identical
GOE spectra at each stationary point, the approach introduced here is overkill.
In Section \ref{sec:multispherical} we apply the methods to a multispherical
spin glass, which is still Gaussian but has a non-GOE spectrum that can vary
between stationary points. Finally, in Section \ref{sec:least.squares} we analyze a model of the sum of squares of random functions, which is non-Gaussian and whose Hessian statistics depend on the conditioning of the energy and gradient.

\subsection{Spherical spin glasses}
\label{sec:ex.spherical}

The spherical spin glasses are a family of models that encompass every
isotropic Gaussian field on the hypersphere defined by all $\mathbf x\in\mathbb R^N$ such that $0=\mathbf x^T\mathbf x-N$. One can consider the models as defined by centered Gaussian functions $H$ such that the covariance between two points in the configuration space is
\begin{equation}
  \overline{H(\mathbf x)H(\mathbf x')}=Nf\left(\frac{\mathbf x^T\mathbf x'}N\right)
\end{equation}
for some function $f$ with positive series coefficients. Such functions can be considered to be made up of all-to-all tensorial interactions, with
\begin{equation}
  H(\mathbf x)
  =\sum_{p=0}^\infty\frac1{p!}\sqrt{\frac{f^{(p)}(0)}{N^{p-1}}}
  \sum_{i_1\cdots i_p}^NJ_{i_1\cdots i_p}x_{i_1}\cdots x_{i_p}
\end{equation}
and the elements of the tensors $J$ being independently distributed with the
unit normal distribution.

The marginal optima of these models can be studied without the methods
introduced in this paper, and have been in the past \cite{Folena_2020_Rethinking,
Kent-Dobias_2023_How}. First, these models are Gaussian, so at large $N$ the
Hessian is statistically independent of the gradient and energy
\cite{Bray_2007_Statistics}. Therefore, conditioning the Hessian can be done
mostly independently from the problem of counting stationary points. Second, in
these models the Hessian at every point in the landscape belongs to the GOE
class with the same width of the spectrum $\mu_\mathrm m=2\sqrt{f''(1)}$.
Therefore, all marginal optima in these systems have the same constant shift
$\mu=\pm\mu_\mathrm m$. Despite the fact the complexity of marginal optima is
well known by simpler methods, it is instructive to carry through the
calculation for this case, since we will learn something about its application in
more nontrivial settings.

The procedure to treat the complexity of the spherical models has been made in
detail elsewhere \cite{Kent-Dobias_2023_How}. Here we make only a sketch of the
steps involved. First the substitutions \eqref{eq:delta.grad},
\eqref{eq:delta.energy}, and \eqref{eq:delta.eigen} are made to convert the
Dirac $\delta$ functions into exponential integrals, and the substitution
\eqref{eq:determinant} is made to likewise convert the determinant.

Once these substitutions have been made, the entire expression
\eqref{eq:min.complexity.expanded} is an exponential integral whose argument is
a linear functional of $H$. This allows for the average to be taken over the
disorder. If we gather all the $H$-dependant pieces into the linear functional
$\mathcal O$ then the average gives
\begin{equation}
  \begin{aligned}
    \overline{
      e^{\sum_a^n\mathcal O_aH(\mathbf x_a)}
    }
    &=e^{\frac12\sum_a^n\sum_b^n\mathcal O_a\mathcal O_b\overline{H(\mathbf x_a)H(\mathbf x_b)}} \\
    &=e^{N\frac12\sum_a^n\sum_b^n\mathcal O_a\mathcal O_bf\big(\frac{\mathbf x_a^T\mathbf x_b}N\big)}
  \end{aligned}
\end{equation}
The result is an integral that only depends on the many vector variables we
have introduced through their scalar products with each other. We therefore make a change of variables in the integration from those vectors to matrices that encode their possible scalar products. These matrices are
\begin{equation} \label{eq:order.parameters}
  \begin{aligned}
    &C_{ab}=\frac1N\mathbf x_a\cdot\mathbf x_b
    \qquad\qquad
    &R_{ab}=-i\frac1N\mathbf x_a\cdot\hat{\mathbf x}_b&
  \\
  &D_{ab}=\frac1N\hat{\mathbf x}_a\cdot\hat{\mathbf x}_b
  &G_{ab}=\frac1N\bar{\pmb\eta}_a^T\pmb\eta_b&
  \\
  &Q_{ab}^{\alpha\gamma}=\frac1N\mathbf s_a^\alpha\cdot\mathbf s_b^\gamma
  &X^\alpha_{ab}=\frac1N\mathbf x_a\cdot\mathbf s_b^\alpha&
  \\
  &\hat X^\alpha_{ab}=-i\frac1N\hat{\mathbf x}_a\cdot\mathbf s_b^\alpha&&
  \end{aligned}
\end{equation}
Order parameters that mix the normal and Grassmann variables generically vanish
in these settings and we don't consider them here \cite{Kurchan_1992_Supersymmetry}.
This transformation changes the measure of the integral, with
\begin{equation}
  \begin{aligned}
    &\prod_{a=1}^nd\mathbf x_a\,\frac{d\hat{\mathbf x}_a}{(2\pi)^N}\,d\bar{\pmb\eta}_a\,d\pmb\eta\,\prod_{\alpha=1}^{m_a}d\mathbf s_a^\alpha \\
    &\quad=dC\,dR\,dD\,dG\,dQ\,dX\,d\hat X\,(\det J)^{N/2}(\det G)^{-N}
  \end{aligned}
\end{equation}
where $J$ is the Jacobian of the transformation and takes the form
\begin{equation} \label{eq:coordinate.jacobian}
  J=\begin{bmatrix}
    C&iR&X_1&\cdots&X_n \\
    iR&D&i\hat X_1&\cdots&i\hat X_n\\
    X_1^T&i\hat X_1^T&Q_{11}&\cdots&Q_{1n}\\
    \vdots&\vdots&\vdots&\ddots&\vdots\\
    X_n^T&i\hat X_n^T&Q_{n1}&\cdots&Q_{nn}
  \end{bmatrix}
\end{equation}
and the contribution of the Grassmann integrals produces its own inverted
Jacobian. The block matrices indicated above are such that $A_{ab}$ is an
$m_a\times m_b$ matrix indexed by the upper indices, while $X_a$ is an $n\times
m_a$ matrix with one lower and one upper index.

After these steps, which follow identically to those more carefully outlined in
the cited papers \cite{Folena_2020_Rethinking, Kent-Dobias_2023_How}, we arrive at a form of the integral as over an effective action
\begin{widetext}
  \begin{equation} \label{eq:spherical.complexity}
  \begin{aligned}
    &\Sigma_{\lambda^*}(E,\mu)
    =\lim_{\beta\to\infty}\lim_{n\to0}\lim_{m_1\cdots m_n\to0}
    \frac1N\frac\partial{\partial n}
    \int dC\,dR\,dD\,dG\,dQ\,dX\,d\hat X\,d\hat\beta\,d\hat\lambda\,
    \exp\Bigg\{
      nN\mathcal S_\mathrm{SSG}(\hat\beta,C,R,D,G\mid E,\mu) \\
      &\qquad
      +nN\mathcal U_\mathrm{SSG}(\hat\lambda,C,Q,X,\hat X\mid\beta)
      +\frac N2\log\det\left[
        I+\begin{bmatrix}
          Q_{11}&\cdots&Q_{1n}\\
          \vdots&\ddots&\vdots\\
          Q_{n1}&\cdots&Q_{nn}
        \end{bmatrix}^{-1}
        \begin{bmatrix}
          X_1^T&i\hat X_1^T\\
          \vdots&\vdots\\
          X_n^T&i\hat X_n^T
        \end{bmatrix}
        \begin{bmatrix}
          C&iR\\iR&D
        \end{bmatrix}^{-1}
        \begin{bmatrix}
          X_1\cdots X_n\\
          i\hat X_1\cdots i\hat X_n
        \end{bmatrix}
      \right]
    \Bigg\}
  \end{aligned}
\end{equation}
where the matrix $J$ is the Jacobian associated with the change of variables
from the $\mathbf x$, $\hat{\mathbf x}$, and $\mathbf s$, and has the form
The structure of the integrand, with the effective action split between two
terms which only share a dependence on the Lagrange multiplier $\omega$ that
enforces the constraint, is generic to Gaussian problems. This is the
appearance in practice of the fact mentioned before that conditions on the
Hessian do not mostly effect the rest of the complexity problem.

The effective action $\mathcal S_\mathrm{SSG}$ is precisely that for the
ordinary complexity of stationary points, or
\begin{equation}
  \begin{aligned}
    &\mathcal S_\mathrm{SSG}(\hat\beta,C,R,D,G\mid E,\mu)
    =\hat\beta E-(r_d+g_d)\mu \\
    &+\frac1n\left\{\frac12\sum_{ab}\left(
      \hat\beta^2f(C_{ab})
      +\big(2\hat\beta R_{ab}-D_{ab}\big)f'(C_{ab})
      +(R_{ab}^2-G_{ab}^2)f''(C_{ab})
    \right)
    +\frac12\log\det\begin{bmatrix}C&iR\\iR^T&D\end{bmatrix}
    -\log\det G\right\}
  \end{aligned}
\end{equation}
where $r_d$ and $g_d$ are the diagonal elements of $R$ and $G$, respectively.
\begin{equation}
  \begin{aligned}
    &\mathcal U_\mathrm{SSG}(\hat\lambda,Q,X,\hat X\mid\lambda^*,\mu,C)
    =\hat\lambda\lambda^*
    +\frac1n\Bigg\{
    \frac12\log\det Q+
      \sum_{a=1}^n\bigg(
        \sum_{\alpha=1}^{m_a}\beta\mu Q_{aa}^{\alpha\alpha}
        +\hat\lambda\mu Q_{aa}^{11}
      \bigg)
    +2\sum_{ab}^nf''(C_{ab})
    \\
     &\qquad\times\Bigg[\beta\sum_\alpha^{m_a}\left(
        \sum_\gamma^{m_b}(Q_{ab}^{\alpha\gamma})^2
        -\hat\beta(X_{ab}^\alpha)^2
        -2X_{ab}^\alpha\hat X_{ab}^\alpha
      \right)
    +\hat\lambda\left(
      \hat\lambda(Q_{ab}^{11})^2
      -\hat\beta(X_{ab}^1)^2
      -2X_{ab}^1\hat X_{ab}^1
    \right)
    +\beta\hat\lambda\left(
        \sum_\alpha^{m_a} Q_{ab}^{\alpha1}
        +\sum_\alpha^{m_b} Q_{ab}^{1\alpha}
    \right)\Bigg]
    \Bigg\}
  \end{aligned}
\end{equation}
\end{widetext}
There are some dramatic simplifications that emerge from the structure of this
particular problem. First, notice that (outside of the `volume' term due to
$J$) the dependence on the parameters $X$ and $\hat X$ are purely quadratic.
Therefore, there will always be a saddle point condition where they are both
zero. In this case, we except this solution to be correct. We can reason about
why this is so: $X$, for instance, quantifies the correlation between the
typical position of stationary points and the direction of their typical
eigenvectors. In an isotropic landscape, where no direction is any more
important than any other, we don't expect such correlations to be nonzero:
where a state is location does not give any information as to the orientation
of its soft directions. On the other hand, in the spiked case, or with an
external field, the preferred direction can polarize both the direction of
typical stationary points \emph{and} their soft eigenvectors. Therefore, in
these instances one must account for solutions with nonzero $X$ and $\hat X$.



When we take $X=\hat X=0$, $Q^{\alpha\beta}_{ab}=\delta_{ab}Q^{\alpha\beta}$
independent, and $Q$ to have the planted replica symmetric form of
\eqref{eq:Q.structure}, we find that
\begin{equation}
  \mathcal U_\mathrm{SSG}(\hat\lambda,Q,0,0\mid\beta,\lambda^*,\mu,C)
  =\mathcal U_\mathrm{GOE}(\hat\lambda,Q\mid\mu,\lambda^*,\beta)
\end{equation}
with $\sigma^2=f''(1)$. That is, the effective action for the terms related to
fixing the eigenvalue in the spherical Kac--Rice problem is exactly the same as
that for the \textrm{GOE} problem. This is perhaps not so surprising, since we
established from the beginning that the Hession of the spherical spin glasses
belongs to the GOE class.

\begin{equation}
  \Sigma_{\lambda^*}(E,\mu)
  =\Sigma(E,\mu)+G_{\lambda^*}(\mu)
\end{equation}
where $G$ is precisely the function \eqref{eq:goe.large.dev} we found in the
case of a GOE matrix added to an identity, with $\sigma=\sqrt{f''(1)}$. We find the marginal complexity by solving
\begin{equation}
  0
  =\frac\partial{\partial\lambda^*}\Sigma_{\lambda^*}(E,\mu_\mathrm m(E))\bigg|_{\lambda^*=0}
  =\frac\partial{\partial\lambda^*}G_{\lambda^*}(\mu_\mathrm m(E))\bigg|_{\lambda^*=0}
\end{equation}
which gives $\mu_m(E)=2\sqrt{f''(1)}$ independent of $E$, as we presaged above. Since $G_0(\mu_\mathrm m)=0$, this gives finally
\begin{equation}
  \Sigma_\mathrm m(E)
  =\Sigma_0(E,\mu_\mathrm m(E))
  =\Sigma(E,\mu_\mathrm m)
\end{equation}
that the marginal complexity in these models is simply the ordinary complexity
evaluated at a fixed trace of the Hessian.

\subsection{Multispherical spin glasses}
\label{sec:multispherical}

The multispherical models are a simple extension of the spherical ones, where
the configuration space is taken to be the union of more than one hypersphere.
Here we consider the specific case where the configuration space is the union
of two $(N-1)$-spheres, with $\Omega=S^{N-1}\times S^{N-1}$, and where the
energy is given by
\begin{equation}
  H(\mathbf x)=H_1(\mathbf x^{(1)})+H_2(\mathbf x^{(2)})-\epsilon\mathbf x^{(1)}\cdot\mathbf x^{(2)}
\end{equation}
for $\mathbf x=[\mathbf x^{(1)},\mathbf x^{(2)}]$ for components $\mathbf
x^{(1)},\mathbf x^{(2)}\in\mathbb R^N$. Each individual sphere energy $H_s$ is
taken to be a centered Gaussian random function with a covariance given in the
usual spherical way by
\begin{equation}
  \overline{H_i(\pmb\sigma_1)H_j(\pmb\sigma_2)}
  =N\delta_{ij}f_i\left(\frac{\pmb\sigma_1\cdot\pmb\sigma_2}N\right)
\end{equation}
with the functions $f_1$ and $f_2$ not necessarily the same. In this problem,
there is an energetic competition between the independent spin glass energies
on each sphere and their tendency to align or anti-align through the
interaction term.

These models have more often been studied with random fully connected couplings
between the spheres, for which it is possible to also use configuration spaces
involving spheres of different sizes \cite{Subag_2021_TAP, Subag_2023_TAP,
Bates_2022_Crisanti-Sommers, Bates_2022_Free, Huang_2023_Strong,
Huang_2023_Algorithmic, Huang_2024_Optimization}.

Because the energy is Gaussian, properties of the Hessian are once again
statistically independent of those of the energy and gradient. However, unlike
the previous example of the spherical models, the spectrum of the Hessian at
different points in the configuration space has different shapes. This appears
in this problem through the presence of a configuration space defined by
multiple constraints, and therefore multiple Lagrange multipliers are necessary
to ensure they are all fixed. The resulting Lagrangian, gradient, and Hessian are
\begin{align}
  L(\mathbf x)&=H(\mathbf x)
  +\frac12\omega^{(1)}\big(\|\mathbf x^{(1)}\|^2-N\big) \\
              &\qquad\qquad\qquad+\frac12\omega^{(2)}\big(\|\mathbf x^{(2)}\|^2-N\big)
  \\
  \nabla H(\mathbf x,\pmb\omega)
              &=\partial H(\mathbf x)+\begin{bmatrix}
    \omega^{(1)}\mathbf x^{(1)} \\
    \omega^{(2)}\mathbf x^{(2)}
  \end{bmatrix}
  \\
  \operatorname{Hess}H(\mathbf x,\pmb\omega)
              &=\partial\partial H(\mathbf x)+\begin{bmatrix}
    \omega^{(1)}I&0 \\
    0&\omega^{(2)}I
  \end{bmatrix}
\end{align}
Like in the spherical model, fixing the trace of the Hessian to $\mu$ is
equivalent to a constraint on the Lagrange multipliers. However, in this case
it corresponds to $\mu=\omega^{(1)}+\omega^{(2)}$, and therefore they are not
uniquely fixed by the trace.

Since the energy in the multispherical models is Gaussian, the properties of
the matrix $\partial\partial H$ are again independent of the energy and
gradient. This means that the form of the Hessian is parameterized solely by
the values of the Lagrange multipliers $\omega^{(1)}$ and $\omega^{(2)}$, just
as $\mu=\omega$ alone parameterized the Hessian in the spherical spin glasses.
Unlike that case, however, the Hessian takes different shapes with different
spectral widths depending on their precise combination. In
Appendix~\ref{sec:multispherical.spectrum} we derive a variational form for the
spectral density of the Hessian in these models using standard methods.

Because of the independence of the Hessian, the method introduced in this
article is not necessary to characterize the marginal minima of this system.
Rather, we could take the spectral density derived in
Appendix~\ref{sec:multispherical.spectrum} and found the Lagrange multipliers
$\omega_1$ and $\omega_2$ corresponding with marginality by tuning the edge of
the spectrum to zero. In some ways the current method is more convenient than
this, since it is a purely variational method and therefore can be reduced to a
since root-finding exercise.

Unlike the constraints on the configurations $\mathbf x$, the constraint on the
tangent vectors $\mathbf s=[\mathbf s^{(1)},\mathbf s^{(2)}]\in\mathbb R^{2N}$
remains the same spherical constraint as before, which implies $N=\|\mathbf
s\|^2=\|\mathbf s^{(1)}\|^2+\|\mathbf s^{(2)}\|^2$. Defining intra- and inter-sphere overlap matrices
\begin{equation}
  Q^{ij,\alpha\gamma}_{ab}=\frac1N\mathbf s^{(i),\alpha}_a\cdot\mathbf s^{(j),\gamma}_b
\end{equation}
this problem no longer has the property that the diagonal of the $Q$s is one,
but instead that $1=Q^{11,\alpha\alpha}_{aa}+Q^{22,\alpha\alpha}_{aa}$. This is
the manifestation of the fact that a normalized vector in the tangent space of the
multispherical model need not be equally spread on the two subspaces, but can
be concentrated in one or the other.

The calculation of the marginal complexity in this problem follows very closely
to that of the spherical spin glasses in the previous subsection, making
immediately the simplifying assumptions that the soft directions of different
stationary points are typically uncorrelated and therefore $X=\hat X=0$ and the
overlaps $Q$ between eigenvectors are only nonzero when in the same replica.
The result has the schematic form of \eqref{eq:spherical.complexity}, but with
different effective actions depending now on overlaps inside each of the two
spheres and between the two spheres. These are
\begin{widetext}
\begin{equation}
  \begin{aligned}
    &\mathcal S_\mathrm{MSG}(\hat\beta,C^{11},R^{11},D^{11},G^{11},C^{22},R^{22},D^{22},G^{22},C^{12},R^{12},R^{21},D^{12},G^{12}
    \mid E,\omega_1,\omega_2)= \\
    &\quad
    \mathcal S_\mathrm{SSG}(\hat\beta,C^{11},R^{11},D^{11},G^{11}\mid E_1,\omega_1)
    +\mathcal S_\mathrm{SSG}(\hat\beta,C^{22},R^{22},D^{22},G^{22}\mid E_2,\omega_2)
    -\epsilon(r^{12}_d+r^{21}_d)+\hat\beta(E-E_1-E_2-\epsilon c_d^{12}) \\
    &\quad
    +\frac12\log\det\left(
    I+
    \begin{bmatrix}C^{11}&iR^{11}\\iR^{11}&D^{11}\end{bmatrix}^{-1}
    \begin{bmatrix}
      C^{12} & iR^{12} \\ iR^{21} & D^{12}
    \end{bmatrix}
    \begin{bmatrix}C^{22}&iR^{22}\\iR^{22}&D^{22}\end{bmatrix}^{-1}
    \begin{bmatrix}
      C^{12} & iR^{21} \\ iR^{21} & D^{12}
    \end{bmatrix}
    \right)
    -\log\det(I+(G^{11}G^{22})^{-1}G^{12}G^{21})
  \end{aligned}
\end{equation}
and
\begin{equation}
  \begin{aligned}
    &\mathcal U_\mathrm{MSG}(\hat q,\hat\lambda,Q^{11},Q^{22},Q^{12}\mid\lambda^*,\omega_1,\omega_2,\beta) \\
    &\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]
    +\hat\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\hat\lambda\sum_a^n(Q^{ii}_{1a})^2+\hat\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}
where again the problem of fixing marginality has completely separated from
that of the complexity. The biggest change between this problem and the
spherical one is that now the spherical constraint in the tangent space at each
stationary point gives the constraint on the order parameters
$q^{11}_d+q^{22}_d=1$. Therefore, the diagonal of the $Q$ matrices cannot be
taken to be 1 as before. To solve the marginal problem, we take each of the
matrices $Q^{11}$, $Q^{22}$, and $Q^{12}$ to have the planted replica symmetric
form \eqref{eq:Q.structure}, but with the diagonal not necessarily equal to 1, so
\begin{equation}
 Q^{ij}=\begin{bmatrix}
   \tilde q^{ij}_d & \tilde q^{ij}_0 & \tilde q^{ij}_0 & \cdots & \tilde q^{ij}_0 \\
   \tilde q^{ij}_0 & q^{ij}_d & q^{ij}_0 & \cdots & q^{ij}_0 \\
   \tilde q^{ij}_0 & q^{ij}_0 & q^{ij}_d & \ddots & q^{ij}_0 \\
   \vdots & \vdots & \ddots & \ddots & \vdots \\
   \tilde q^{ij}_0 & q^{ij}_0 & q^{ij}_0 & \cdots & q^{ij}_d
 \end{bmatrix}
\end{equation}

\begin{widetext}
  \begin{equation}
    \begin{aligned}
      &\sum_{i=1,2}f_i''(1)\left[
        \beta^2\left(
          (\tilde q^{ii}_d)^2
          -(q^{ii}_d)^2
          +2(q^{ii}_0)^2
          -2(\tilde q^{ii}_0)^2
        \right)
        -2\beta\hat\lambda\left(
          (\tilde q^{ii}_d)^2-(\tilde q^{ii}_0))^2
        \right)
        +\hat\lambda^2(\tilde q^{ii}_d)^2
      \right]
      +\hat\lambda\left(
        \tilde q^{11}_d\omega_1+\tilde q^{22}_d\omega_2+2\tilde q^{12}_d
      \right) \\
      &-\beta\left(
        (\tilde q^{11}_d-q^{11}_d)\omega_1
        +(\tilde q^{22}_d-q^{22}_d)\omega_2
        -2\epsilon(\tilde q^{12}_d-q^{12}_d)
      \right) \\
      &+\frac12\log\bigg[
        \left(
          2q^{12}_0\tilde q^{12}_0-\tilde q^{12}_0(\tilde q^{12}_d+q^{12}_d)
          -2\tilde q^{11}_0q^{22}_0+\tilde q^{11}_d\tilde q^{22}_0+\tilde q^{11}_0q^{22}_d
        \right)
        \left(
          2q^{12}_0\tilde q^{12}_0-\tilde q^{12}_0(\tilde q^{12}_d+q^{12}_d)
          -2q^{11}_0\tilde q^{22}_0+q^{11}_d\tilde q^{22}_0+\tilde q^{11}_0\tilde q^{22}_d
        \right) \\
      &\qquad\qquad+2\left(3(q^{12}_0)^2-(\tilde q^{12}_0)^2-2q^{12}_0q^{12}_d-3q^{11}_0q^{22}_0+q^{11}_dq^{22}_0+\tilde q^{11}_0\tilde q^{22}_0+q^{11}_0q^{22}_d
        \right)\left(
          (\tilde q^{12}_0)^2-(\tilde q^{12}_d)^2-\tilde q^{11}_0\tilde q^{22}_0+\tilde q^{11}_d\tilde q^{22}_d
        \right) \\
      &\qquad\qquad+\left(
        2(q^{12}_0)^2-(\tilde q^{12}_0)^2-(q^{12}_d)^2-2q^{11}_0q^{22}_0+\tilde q^{11}_0\tilde q^{22}_0+q^{11}_dq^{22}_d
      \right)\left(
        (\tilde q^{12}_0)^2-(\tilde q^{12}_d)^2-\tilde q^{11}_0\tilde q^{22}_0+\tilde q^{11}_d\tilde q^{22}_d
      \right)
      \bigg]
      \\
      &-\log\left[(q^{11}_d-q^{11}_0)(q^{22}_d-q^{22}_0)-(q^{12}_d-q^{12}_0)^2\right]
    \end{aligned}
  \end{equation}
\end{widetext}
To make the limit to zero temperature, we once again need an ansatz for the
asymptotic behavior of the overlaps. These take the form
$q^{ij}_0=q^{ij}_d-y^{ij}_0\beta^{-1}-z^{ij}_0\beta^{-2}$, with the same for
the tilde variables. Notice that in this case, the asymptotic behavior of the
off diagonal elements is to approach the value of the diagonal rather than one.
We also require $\tilde q^{ij}_d=q^{ij}_d-\tilde y^{ij}_d\beta^{-1}-\tilde
z^{ij}_d\beta^{-2}$, i.e., that the tilde diagonal term also approaches the
same diagonal value.

As before, in order for the volume term to stay finite, there are necessary
constraints on the values $y$. These are
\begin{align}
  \frac12(y^{11}_d-\tilde y^{11}_d)=y^{11}_0-\tilde y^{11}_0 \\
  \frac12(y^{22}_d-\tilde y^{22}_d)=y^{22}_0-\tilde y^{22}_0 \\
  \frac12(y^{12}_d-\tilde y^{12}_d)=y^{12}_0-\tilde y^{12}_0
\end{align}
One can see that when the diagonal elements are all equal, this requires the
$y$s for the off-diagonal elements to be equal, as in the GOE case. Here, since
the diagonal elements are not necessarily equal, we have a more general
relationship.

\begin{figure}
  \includegraphics{figs/msg_marg_legend.pdf}

  \includegraphics{figs/msg_marg_params.pdf}
  \hfill
  \includegraphics{figs/msg_marg_spectra.pdf}

  \caption{
    \textsc{Left}: Values of the Lagrange multipliers $\omega_1$ and $\omega_2$
    corresponding to a marginal spectrum for multispherical spin glasses with
    $\sigma_1^2=f_1''(1)=1$, $\sigma_2^2=f_2''(1)=1$, and various $\epsilon$.
    \textsc{Right}: Spectra corresponding to the parameters $\omega_1$ and
    $\omega_2$ marked by the circles on the lefthand plot.
  } \label{fig:msg.marg}
\end{figure}

Fig.~\ref{fig:msg.marg} shows the examples of the Lagrange multipliers
necessary for marginality in a set of multispherical spin glasses at various
couplings $\epsilon$, along with some of the corresponding spectra. As
expected, the method correctly picks out values of the Lagrange multipliers
that result in marginal spectra.

Multispherical spin glasses may be an interesting platform for testing ideas
about which among the possible marginal minima actually attract the dynamics,
and which do not. In the limit where $\epsilon=0$ and the configurations of the
two spheres are independent, the minima found should be marginal on both
sphere's energies. Just because technically on the expanded configuration space
a deep and stable minimum on one sphere and a marginal minimum on the other is
a marginal minimum on the whole space doesn't mean the deep and stable minimum
is any easier to find. This intuitive idea that is precise in the zero-coupling
limit should continue to hold at small nonzero coupling, and perhaps reveal
something about the inherent properties of marginal minima that do not tend to be found
by algorithms.

\subsection{Random nonlinear least squares}
\label{sec: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 minus 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.

\cite{Urbani_2023_A, Kamali_2023_Dynamical, Kamali_2023_Stochastic, Urbani_2024_Statistical}
\cite{Montanari_2023_Solving, Montanari_2024_On}
\cite{Subag_2020_Following}

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}
As in the spherical and multispherical models, fixing the trace of the Hessian
at largest order in $N$ is equivalent to constraining the value of the Lagrange
multiplier $\omega=\mu$, since the trace of the random parts of the Hessian
matrix contribute typical values at a lower order in $N$.

The derivation of the marginal complexity for this model is complicated, but
can be made schematically like that of the derivation of the equilibrium free
energy by use of superspace coordinates \cite{DeWitt_1992_Supermanifolds}.
The use of superspace coordinates in the geometry and dynamics of disordered
systems is well-established. Here, we introduce a novel extension of the
traditional approach to incorporate the marginality condition.
Consider supervectors in the $\mathbb R^{N|4}$ superspace of the form
\begin{equation}
  \pmb\phi_a^\alpha(1,2)
  =\mathbf x_a
  +\bar\theta_1\pmb\eta_a+\bar{\pmb\eta}_a\theta_1
  +i\hat{\mathbf x}_a\bar\theta_1\theta_1
  +\mathbf s_a^\alpha(\bar\theta_1\theta_2+\bar\theta_2\theta_1)
\end{equation}
The traditional complexity problem, outlined in the appendix
\ref{sec:dominant.complexity}, involves a supervector without the last term.
\begin{widetext}
  Using the notation of Section \ref{sec:superspace_kac-rice}, the replicated
  number of stationary points conditioned on energy $E$, trace $\mu$, and
  minimum eigenvalue $\lambda^*$ is then given by
\begin{equation}
  \begin{aligned}
    \mathcal N(E,\mu,\lambda^*)^n
    &=\int\prod_{a=1}^n\lim_{m_a\to0}\prod_{\alpha=1}^{m_a}d\pmb\phi_a^\alpha
    \exp\left\{
      \delta^{\alpha1}N(\hat\beta_aE+\hat\lambda_a\lambda^*)
      -\frac12\int d1\,d2\,B_a^\alpha(1,2)\left[\sum_{k=1}^MV_k(\pmb\phi_a^\alpha)^2
      -\mu(\|\pmb\phi_a^\alpha\|^2-N)\right]
    \right\}
  \end{aligned}
\end{equation}
The first step to evaluate this expression is to linearize the dependence on the random functions $V$. This is accomplished by inserting into the integral a Dirac $\delta$ function fixing the value of the energy for each replica, or
\begin{equation}
  \delta\big(
    V_k(\pmb\phi_a^\alpha(1,2))-v_{ka}^\alpha(1,2)
  \big)
  =
  \int d\hat v_{ka}^\alpha\exp\left[
    i\int d1\,d2\,\hat v_{ka}^\alpha(1,2)
    \big(V_k(\pmb\phi_a^\alpha(1,2))-v_{ka}^\alpha(1,2)\big)
  \right]
\end{equation}
where we have introduced auxiliary fields $\hat v$. With this inserted into the
integral, all other instances of $V$ are replaced by $v$, and the only
remaining dependence on the disorder is from the term $\hat vV$ arising from
the Fourier representation of the Dirac $\delta$ function. This term is linear in $V$, and therefore the random functions can be averaged over to produce
\begin{equation}
  \overline{
    \exp\left[
      i\sum_{ka\alpha}\int d1\,d2\,\hat v_{ka}^\alpha(1,2)
      V_k(\pmb\phi_a^\alpha(1,2))
    \right]
  }
  =
  -\frac N2\sum_{ab}^n\sum_{\alpha\gamma}^{m_a}\sum_k^{\alpha N}\int d1\,d2\,d3\,d4\,
  \hat v_{ka}^\alpha(1,2)f\big(\pmb\phi_a^\alpha(1,2)^T\pmb\phi_b^\gamma(3,4)\big)\hat v_{kb}^\gamma(3,4)
\end{equation}
The entire integrand is now quadratic in the $v$ and $\hat v$ with the kernel
\begin{equation}
  \begin{bmatrix}
    B_a^\alpha(1,2)\delta(1,3)\delta(2,4)\delta_{ab}\delta^{\alpha\gamma} & i\delta(1,3)\,\delta(2,4) \delta_{ab}\delta^{\alpha\gamma}\\
    i\delta(1,3)\,\delta(2,4) \delta_{ab}\delta^{\alpha\gamma}& f\big(\pmb\phi_a^\alpha(1,2)^T\pmb\phi_b^\gamma(3,4)\big)
  \end{bmatrix}
\end{equation}
The integration over the $v$ and $\hat v$ results in a term in the effective action of the form
\begin{equation}
  -\frac M2\log\operatorname{sdet}\left(
    \delta(1,3)\,\delta(2,4) \delta_{ab}\delta^{\alpha\gamma}
    +B_a^\alpha(1,2)f\big(\pmb\phi_a^\alpha(1,2)^T\pmb\phi_b^\gamma(3,4)\big)
  \right)
\end{equation}
When expanded, this supermatrix is constructed of the scalar products of the
real and Grassmann vectors that make up $\pmb\phi$. The change of variables to
these order parameters again results in the Jacobian of \eqref{eq:coordinate.jacobian}, contributing
\begin{equation}
  \frac N2\log\det J(C,R,D,Q,X,\hat X)-\frac N2\log\det G^2
\end{equation}
Up to this point, the expressions above are general and independent of a given
ansatz. However, we expect that the order parameters $X$ and $\hat X$ are zero,
since this case is isotropic. Applying this ansatz here avoids a dramatically
more complicated expression for the effective action found in the case with
arbitrary $X$ and $\hat X$. We also will apply the ansatz that $Q_{ab}^{\alpha\gamma}$ is zero for $a\neq b$, which is equivalent to assuming that the soft
directions of typical pairs of stationary points are uncorrelated, and further
that $Q^{\alpha\gamma}=Q_{aa}^{\alpha\gamma}$ independently of the index $a$,
implying that correlations in the tangent space of typical stationary points
are the same.

Given these simplifying forms of the ansatz, taking the superdeterminant yields
\begin{equation}
  \begin{aligned}
    \log\det\left\{
      \left[
        f'(C)\odot D-\hat\beta I+\left(R^{\circ2}-G^{\circ2}+I\sum_{\alpha\gamma}2(\delta^{\alpha1}\hat\lambda+\beta)(\delta^{\gamma1}\hat\lambda+\beta)Q_{\alpha\gamma}^2\right)\odot f''(C)
      \right]f(C)
      +(I-R\odot f'(C))^2
    \right\} \\
    +n\log\det_{\alpha\gamma}(\delta_{\alpha\gamma}-2(\delta_{\alpha1}\hat\lambda+\beta)Q_{\alpha\gamma})
    -2\log\det(I+G\odot f'(C))
  \end{aligned}
\end{equation}
where once again $\odot$ is the Hadamard product and $A^{\circ n}$ gives the
Hadamard power of $A$. We can already see one substantive difference between
the structure of this problem and that of the spherical models: the effective
action in this case mixes the order parameters $G$ due to the Grassmann variables with the
ones $C$, $R$, and $D$ due to the other variables. Notice further that the dependence on $Q$ due to the marginal constraint is likewise no longer separable. This is the realization of
the fact that the Hessian properties are no longer independent of the energy
and gradient.

Now we have reduced the problem to an extremal one over the order parameters
$\hat\beta$, $\hat\lambda$, $C$, $R$, $D$, $G$, and $Q$, it is time to make an
ansatz for the form of order we expect to find. We will focus on a regime where
the structure of stationary points is replica symmetric, and further where
typical pairs of stationary points have no overlap. This requires that $f(0)=0$, or that there is no constant term in the random polynomials. This gives
\begin{align}
  C=I && R=rI && D = dI && G = gI
\end{align}
We further take a planted replica symmetric structure for the matrix $Q$,
identical to that in \eqref{eq:Q.structure}. The resulting effective action is
the same as if we had made an annealed calculation in the complexity, though
the previous expressions are general.
\begin{equation}
  \begin{aligned}
    \mathcal S_\beta
    =\hat\beta E-\mu(r+g)
    +\frac12\log\frac{d+r^2}{g^2}\frac{1-2q_0+\tilde q_0^2}{(1-q_0)^2}
    -\frac\alpha2\log\left(\frac{1-f'(1)(2\beta(1-q_0)+\hat\lambda-(1-2q_0+\tilde q_0^2)\beta(\beta+\hat\lambda)f'(1))}{(1-(1-q_0)\beta f'(1))^2}\right)
    \\
    -\frac12\mu\hat\lambda+\hat\lambda\lambda^*-\frac\alpha2\log\left[
      \frac{
        \big[f'(1)d-\hat\beta-f''(1)(r^2-g^2+q_0^2\beta^2-\tilde q_0^2\beta(\beta+\hat\lambda)+\beta\hat\lambda+\frac12\hat\lambda^2)\big]f(1)+(1-rf'(1))^2
      }{
        (1+gf'(1))^2
      }
    \right]
  \end{aligned}
\end{equation}
We expect as before the limits of $q_0$ and $\tilde q_0$ as $\beta$ goes to
infinity to approach one, defining their asymptotic expansion as in
\eqref{eq:q0.limit} and \eqref{eq:q0t.limit}. Upon making this substitution and
taking the zero-temperature limit, we find
\begin{equation}
  \begin{aligned}
    \mathcal S_\infty
    =\hat\beta E-\mu(r+g)
    +\frac12\log\frac{d+r^2}{g^2}\frac{y_0^2-\Delta z}{y_0^2}
    -\frac\alpha2\log\left(
      \frac{
        1-(2y_0+\hat\lambda)f'(1)+(y_0^2-\Delta z)f'(1)^2
      }{(1-y_0f'(1))^2}
    \right)
    \\
    -\frac12\mu\hat\lambda+\hat\lambda\lambda^*-\frac\alpha2\log\left[
      \frac{
        \big[f'(1)d-\hat\beta-f''(1)(r^2-g^2+2y_0\hat\lambda+\Delta z+\frac12\hat\lambda^2)\big]f(1)+\big[1-rf'(1)\big]^2
      }{
        (1+gf'(1))^2
      }
    \right]
  \end{aligned}
\end{equation}
\begin{equation}
  \Sigma_{\lambda^*}(E,\mu)=\operatorname{extremum}_{\hat\beta,r,d,g,y_0,\Delta z,\hat\lambda}\mathcal S_\infty
\end{equation}
\end{widetext}

\begin{figure}
  \includegraphics{figs/most_squares_complexity.pdf}
  \caption{
    Dominant and marginal complexity in the nonlinear sum of squares problem
    for $\alpha=\frac32$ and $f(q)=q^2+q^3$. The ground state energy
    $E_\mathrm{gs}$ and the threshold energy $E_\mathrm{th}$ are marked on the
    plot.
  }
\end{figure}

\section{Conclusions}
\label{sec:conclusion}

We have introduced a method for conditioning complexity on the marginality of
stationary points. This method is in principal completely general, and permits
this conditioning without first needing to understand the entire Hessian
statistics. We used our approach to study the marginal complexity in three
different models of random landscapes, showing that the method works and can be
applied to models whose marginal complexity was not previously known. In our
companion paper, we further show that the marginal complexity in the third
model of random nonlinear least squares can be used to effectively bound
algorithmic performance \cite{Kent-Dobias_2024_Algorithm-independent}.

There are some limitations to the approach we have largely relied in this
paper. The main limitation is our restriction to signalless landscapes, where
there is no symmetry-breaking favored direction. This allowed us to neglect the
presence of stationary points with isolated eigenvalues as atypical, and
therefore apply the marginal conditioning using a variational principle.
However, most models of interest in inference have a nonzero signal strength
and therefore often have typical stationary points with an isolated eigenvalue.
As we described earlier, marginal complexity can still be analyzed in these
systems by tuning the shift $\mu$ until the large-deviation principle breaks
down and an imaginary part of the complexity appears. However, this is an
inconvenient measure. It's possible that a variational approach can be
preserved by treating the direction toward and the directions orthogonal to the
signal differently. This problem merits further research.

\begin{acknowledgements}
  JK-D is supported by a \textsc{DynSysMath} Specific Initiative of the INFN.
\end{acknowledgements}

\appendix

\section{A primer on superspace}
\label{sec:superspace}

The superspace $\mathbb R^{N|2D}$ is a vector space with $N$ real indices and
$2D$ Grassmann indices $\bar\theta_1,\theta_1,\ldots,\bar\theta_D,\theta_D$.
The Grassmann indices anticommute like fermions. Their integration is defined by
\begin{equation}
  \int d\theta\,\theta=1
  \qquad
  \int d\theta\,1=0
\end{equation}
Because the Grassmann indices anticommute, their square is always zero.
Therefore, any series expansion of a function with respect to a given Grassmann
index will terminate exactly at linear order, while a series expansion with
respect to $n$ Grassmann variables will terminate exactly at $n$th order. If
$f$ is an arbitrary function, then
\begin{equation}
  \int d\theta\,f(a+b\theta)
  =\int d\theta\,\left[f(a)+f'(a)b\theta\right]
  =f'(a)b
\end{equation}
This kind of behavior of integrals over the Grassmann indices makes them useful
for compactly expressing the Kac--Rice measure. To see why, consider the
specific superspace $\mathbb R^{N|2}$, where an arbitrary vector can be expressed as
\begin{equation}
  \pmb\phi(1)=\mathbf x+\bar\theta_1\pmb\eta+\bar{\pmb\eta}\theta_1+\bar\theta_1\theta_1i\hat{\mathbf x}
\end{equation}
where $\mathbf x,\hat{\mathbf x}\in\mathbb R^N$ and $\bar{\pmb\eta},\pmb\eta$ are
$N$-dimensional Grassmann vectors. The dependence of $\pmb\phi$ on 1 indicates
the index of Grassmann variables $\bar\theta_1,\theta_1$ inside, since we will
sometimes want to use, e.g., $\pmb\phi(2)$ defined identically save for
substitution by $\bar\theta_2,\theta_2$. Consider the series expansion of an arbitrary function $f$ of this supervector:
\begin{equation}
  \begin{aligned}
    f\big(\pmb\phi(1)\big)
    &=f(\mathbf x)
    +\big(\bar\theta_1\pmb\eta+\bar{\pmb\eta}\theta_1+\bar\theta_1\theta_1i\hat{\mathbf x}\big)^T\partial f(\mathbf x) \\
    &\quad+\frac12\big(\bar\theta_1\pmb\eta+\bar{\pmb\eta}\theta_1\big)^T\partial\partial f(\mathbf x)\big(\bar\theta_1\pmb\eta+\bar{\pmb\eta}\theta_1\big) \\
    &=f(\mathbf x)
    +\big(\bar\theta_1\pmb\eta+\bar{\pmb\eta}\theta_1+\bar\theta_1\theta_1i\hat{\mathbf x}\big)^T\partial f(\mathbf x) \\
    &\qquad-\bar\theta_1\theta_1\bar{\pmb\eta}^T\partial\partial f(\mathbf x)\pmb\eta
  \end{aligned}
\end{equation}
where the last step we used the fact that the Hessian matrix is symmetric and
that squares of Grassmann indicies vanish. Using the integration rules defined above, we find
\begin{equation}
  \int d\theta_1\,d\bar\theta_1\,f\big(\pmb\phi(1)\big)
  =i\hat{\mathbf x}^T\partial f(\mathbf x)-\bar{\pmb\eta}^T\partial\partial f(\mathbf x)\pmb\eta
\end{equation}
These two terms are precisely the exponential representation of the Dirac
$\delta$ function of the gradient and determinant of the Hessian (without
absolute value sign) that make up the basic Kac--Rice measure, so that we can write
\begin{equation}
  \begin{aligned}
    &\int d\mathbf x\,\delta\big(\nabla H(\mathbf x)\big)\,\det\operatorname{Hess}H(\mathbf x) \\
    &\qquad=\int d\mathbf x\,d\bar{\pmb\eta}\,d\pmb\eta\,\frac{d\hat{\mathbf x}}{(2\pi)^N}\,e^{i\hat{\mathbf x}^T\nabla H(\mathbf x)-\bar{\pmb\eta}^T\operatorname{Hess}H(\mathbf x)\pmb\eta} \\
    &\qquad=\int d\pmb\phi\,e^{\int d1\,H(\pmb\phi(1))}
  \end{aligned}
\end{equation}
where we have written $d1=d\theta_1\,d\bar\theta_1$ and $d\pmb\phi=d\mathbf
x\,d\bar{\pmb\eta}\,d\pmb\eta\,\frac{d\hat{\mathbf x}}{(2\pi)^N}$. Besides some
deep connections to the physics of BRST, this compact notation dramatically
simplifies the analytical treatment of the problem. The energy of stationary points can also be fixed using this notation, by writing
\begin{equation}
  \int d\pmb\phi\,\frac{d\hat\beta}{2\pi}\,e^{\hat\beta E+\int d1\,(1-\hat\beta\bar\theta_1\theta_1)H(\pmb\phi(1))}
\end{equation}
which a small calculation confirms results in the same expression as \eqref{eq:delta.energy}.

The reason why this simplification is
possible is because there are a large variety of superspace algebraic and
integral operations with direct corollaries to their ordinary real
counterparts. For instance, consider a super linear operator $M(1,2)$, which
like the super vector $\pmb\phi$ is made up of a linear combination of $N\times
N$ regular or Grassmann matrices indexed by every nonvanishing combination of
the Grassmann indices $\bar\theta_1,\theta_1,\bar\theta_2,\theta_2$. Such a supermatrix acts on supervectors by ordinary matrix multiplication and convolution in the Grassmann indices, i.e.,
\begin{equation}
  (M\pmb\phi)(1)=\int d1\,M(1,2)\pmb\phi(2)
\end{equation}
The identity supermatrix is given by
\begin{equation}
  \delta(1,2)=(\bar\theta_1-\bar\theta_2)(\theta_1-\theta_2)I
\end{equation}
Integrals involving superfields contracted into such operators result in schematically familiar expressions, like that of the standard Gaussian:
\begin{equation}
  \int d\pmb\phi\,e^{\int\,d1\,d2\,\pmb\phi(1)^TM(1,2)\pmb\phi(2)}
  =(\operatorname{sdet}M)^{-1/2}
\end{equation}
where the usual role of the determinant is replaced by the superdeterminant.
The superdeterminant can be defined using the ordinary determinant by writing a
block version of the matrix $M$: if $\mathbf e(1)=\{1,\bar\theta_1\theta_1\}$ is
the basis vector of the even subspace of the superspace and $\mathbf
f(1)=\{\bar\theta_1,\theta_1\}$ is that of the odd subspace, then we can form a
block representation of $M$ in analogy to the matrix form of an operator in quantum mechanics by
\begin{equation}
  \int d1\,d2\,\begin{bmatrix}
    \mathbf e(1)M(1,2)\mathbf e(2)^T
    &
    \mathbf e(1)M(1,2)\mathbf f(2)^T
    \\
    \mathbf f(1)M(1,2)\mathbf e(2)^T
    &
    \mathbf f(1)M(1,2)\mathbf f(2)^T
  \end{bmatrix}
  =\begin{bmatrix}
    A & B \\ C & D
  \end{bmatrix}
\end{equation}
where each of the blocks is a $2N\times 2N$ real matrix. Then the
superdeterminant of $M$ is given by
\begin{equation}
  \operatorname{sdet}M=\det(A-BD^{-1}C)\det(D)^{-1}
\end{equation}
which is the same for the normal equation for the determinant of a block matrix
save for the inverse of $\det D$. The same method can be used to calculate the
superdeterminant in arbitrary superspaces, where for $\mathbb R^{N|2D}$ each
basis has $2^{2D-1}$ elements. For instance, for $\mathbb R^{N|4}$ we have $\mathbf e(1,2)=\{1,\bar\theta_1\theta_1,\bar\theta_2\theta_2,\bar\theta_1\theta_2,\bar\theta_2\theta_1,\bar\theta_1\bar\theta_2,\theta_1\theta_2,\bar\theta_1\theta_1\bar\theta_2\theta_2\}$ and $\mathbf f(1,2)=\{\bar\theta_1,\theta_1,\bar\theta_2,\theta_2,\bar\theta_1\theta_1\bar\theta_2,\bar\theta_2\theta_2\theta_1,\bar\theta_1\theta_1\theta_2,\bar\theta_2\theta_2\theta_1\}$.

\section{BRST symmetry}
\label{sec:brst}

The superspace representation is also helpful because it can make manifest an
unusual symmetry in the dominant complexity of minima that would otherwise be
obfuscated. This arises from considering the Kac--Rice formula as a kind of
gauge fixing procedure \cite{Zinn-Justin_2002_Quantum}. Around each stationary
point consider making the coordinate transformation $\mathbf u=\nabla H(\mathbf
x)$. Then in the absence of fixing the trace, the Kac--Rice measure becomes
\begin{equation}
    \int d\nu(\mathbf x,\pmb\omega\mid E)
    =\int\sum_\sigma d\mathbf u\,\delta(\mathbf u)\,
    \delta\big(NE-H(\mathbf x_\sigma)\big)
\end{equation}
where the sum is over stationary points. This integral has a symmetry of its
measure of the form $\mathbf u\mapsto\mathbf u+\delta\mathbf u$. Under the
nonlinear transformation that connects $\mathbf u$ and $\mathbf x$, this
implies a symmetry of the measure in the Kac--Rice integral of $\mathbf
x\mapsto\mathbf x+(\operatorname{Hess}H)^{-1}\delta\mathbf u$. This symmetry, while exact, is
nonlinear and difficult to work with.

When the absolute value sign has been dropped and Grassmann vectors introduced,
this symmetry can be simplified considerably. Due to the expansion properties
of Grassmann integrals, any appearance of $-\bar{\pmb\eta}\pmb\eta^T$ in the
integrand resolves to $(\operatorname{Hess}H)^{-1}$. The
symmetry of the measure can then be written
\begin{equation}
  \mathbf x\mapsto \mathbf x-\bar{\pmb\eta}\pmb\eta^T\delta\mathbf u
  =\mathbf x+\bar{\pmb\eta}\delta\epsilon
\end{equation}
where $\delta\epsilon=-\pmb\eta^T\delta\mathbf u$ is a Grassmann number. This
establishes that $\delta\mathbf x=\bar{\pmb\eta}\delta\epsilon$, now linear. The rest of
the transformation can be built by requiring that the action is invariant after
expansion in $\delta\epsilon$. Ignoring for a moment the piece of the measure
fixing the trace of the Hessian, this gives
\begin{align}
  \delta\mathbf x=\bar{\pmb\eta}\,\delta\epsilon &&
  \delta\hat{\mathbf x}=-i\hat\beta\bar{\pmb\eta}\,\delta\epsilon &&
  \delta\pmb\eta=-i\hat{\mathbf x}\,\delta\epsilon  &&
  \delta\bar{\pmb\eta}=0
\end{align}
so that the differential form of the symmetry is
\begin{equation}
  \mathcal D=\bar{\pmb\eta}\frac\partial{\partial\mathbf x}
  -i\hat\beta\bar{\pmb\eta}\frac\partial{\partial\hat{\mathbf x}}
  -i\hat{\mathbf x}\frac\partial{\partial\pmb\eta}
\end{equation}
The Ward identities associated with this symmetry give rise to relationships among the order parameters. These identities are
\begin{align}
  0=\frac1N\mathcal D\langle\mathbf x_a^T\pmb\eta_b\rangle
  =\frac1N\left[
    \langle\bar{\pmb\eta}_a^T\pmb\eta_b\rangle-
    i\langle\mathbf x_a^T\hat{\mathbf x}_b\rangle
  \right]
  =G_{ab}+R_{ab} \\
  0=\frac iN\mathcal D\langle\hat{\mathbf x}_a^T\pmb\eta_b\rangle
  =\frac1N\left[
    \hat\beta\langle\bar{\pmb\eta}_a^T\pmb\eta_b\rangle
    +\langle\hat{\mathbf x}_a^T\hat{\mathbf x}_b\rangle
  \right]
  =\hat\beta G_{ab}+D_{ab}
\end{align}
These identities establish $G_{ab}=-R_{ab}$ and $D_{ab}=\hat\beta R_{ab}$,
allowing elimination of the matrices $G$ and $D$ in favor of $R$. Fixing the
trace to $\mu$ explicitly breaks this symmetry, and the simplification is lost.

\section{Spectral density in the multispherical spin glass}
\label{sec:multispherical.spectrum}

In this appendix we derive an expression for the asymptotic spectral density in
the two-sphere multispherical spin glass that we describe in Section
\ref{sec:multispherical}. \cite{Livan_2018_Introduction}
\begin{equation}
  \begin{aligned}
    &G(\lambda)
    =\lim_{n\to0}\int\|\mathbf y_1\|^2\,\prod_{a=1}^nd\mathbf y_a\,
    \exp\left\{
      -\frac12\mathbf y_a^T(\operatorname{Hess}H(\mathbf x,\pmb\omega)+\lambda I)\mathbf y_a
    \right\} \\
    &
    =\lim_{n\to0}\int\big(\|\mathbf y_1^{(1)}\|^2+\|\mathbf y_1^{(2)}\|^2\big)\,\prod_{a=1}^nd\mathbf y_a\, \\
    &\times\exp\left\{
      -\frac12\begin{bmatrix}\mathbf y_a^{(1)}\\\mathbf y_a^{(2)}\end{bmatrix}^T
      \left(
        \begin{bmatrix}
          \operatorname{Hess}H_1(\mathbf x^{(1)},\omega_1) & -\epsilon \\
          -\epsilon & \operatorname{Hess}H_2(\mathbf x^{(2)},\omega_2)
        \end{bmatrix}
        +\lambda I
      \right)\begin{bmatrix}\mathbf y_a^{(1)}\\\mathbf y_a^{(2)}\end{bmatrix}
    \right\} \\
  \end{aligned}
\end{equation}
If $Y_{ab}^{(ik)}=\frac1N\mathbf y_a^{(i)}\cdot\mathbf y_b^{(j)}$ is the matrix
of overlaps of the $\mathbf y$, then a short and standard calculation yields
\begin{equation}
  G(\lambda)=N\lim_{n\to0}\int dY\,(Y_{11}^{(11)}+Y_{11}^{(22)})\,
  e^{nN\mathcal S(Y)}
\end{equation}
for
\begin{equation}
  \begin{aligned}
    &\mathcal S(Y)
    =\frac1n\sum_{ab}\left[
      \sigma_1^2(Y_{ab}^{(11)})^2
      +\sigma_2^2(Y_{ab}^{(22)})^2
    \right]+\frac12\log\det\begin{bmatrix}
      Y^{(11)}&Y^{(12)}\\Y^{(12)}&Y^{(22)}
    \end{bmatrix}\\
    &+\frac1n\sum_a^n\left[
      2\epsilon Y_{aa}^{(12)}
      -\omega_1Y_{aa}^{(11)}
      -\omega_2Y_{aa}^{(22)}
      +\lambda(Y_{aa}^{(11)}
      +Y_{aa}^{(22)})
    \right]
  \end{aligned}
\end{equation}
Making the replica symmetric ansatz $Y_{ab}^{(ij)}=y^{(ij)}\delta_{ab}$ yields
\begin{equation}
  \begin{aligned}
    &\mathcal S(y)
    =
      \sigma_1^2(y^{(11)})^2
      +\sigma_2^2(y^{(22)})^2
    +\frac12\log(
      y^{(11)}y^{(22)}-y^{(12)}y^{(12)}
      )\\
    &+2\epsilon y^{(12)}
      -\omega_1y^{(11)}
      -\omega_2y^{(22)}
      +\lambda(y^{(11)}
      +y^{(22)})
  \end{aligned}
\end{equation}
\begin{equation}
  \overline{G(\lambda)}
  =N(y^{(11)}+y^{(22)})
\end{equation}
\begin{equation}
  \rho(\lambda)
  =\frac1{i\pi N}
  \left(
    \overline{G(\lambda+i0^+)}-\overline{G(\lambda+i0^-)}
  \right)
\end{equation}

\section{Complexity of dominant optima in the least-squares problem}
\label{sec:dominant.complexity}

Here we share an outline of the derivation of formulas for the complexity of
dominant optima in the random nonlinear least squares problem of section
\ref{sec:least.squares}. While in this paper we only treat problems with a
replica symmetric structure, formulas for the effective action are generic to
any structure and provide a starting point for analyzing the challenging
full-RSB setting.

Using the $\mathbb R^{N|2}$ superfields
\begin{equation}
  \pmb\phi_a(1)=\mathbf x+\bar\theta_1\pmb\eta+\bar{\pmb\eta}\theta_1+\bar\theta_1\theta_1\hat{\mathbf x},
\end{equation}
the replicated count of stationary points can be written
\begin{equation}
  \begin{aligned}
    &\mathcal N(E,\mu)^n
    =\int\prod_{a=1}^nd\hat\beta_a\,d\pmb\phi_a\,
    \\
    &\qquad\times\exp\left[
      \hat\beta_a E-\frac12\int d1\,B_a(1)\sum_{k=1}^MV^k(\pmb\phi_a(1))^2
    \right]
  \end{aligned}
\end{equation}
for $B_a(1)=1-\hat\beta_a\bar\theta_1\theta_1$.
The derivation of the complexity follows from here nearly identically to that
in Appendix A.2 of \citeauthor{Fyodorov_2022_Optimization} with superoperations
replacing standard ones \cite{Fyodorov_2022_Optimization}. First we insert
Dirac $\delta$ functions to fix each of the $M$ energies $V^k(\pmb\phi_a(1))$ as
\begin{equation} \label{eq:Vv.delta}
  \begin{aligned}
    &\int dv^k_a\,\delta\big(V^k(\pmb\phi_a(1))-v^k_a(1)\big)
    \\
    &\quad=\int dv^k_a\,d\hat v^k_a\,\exp\left[i\int d1\,\hat v^k_a(1)\big(V^k(\pmb\phi_a(1))-v^k_a(1)\big)\right]
  \end{aligned}
\end{equation}
The squared $V^k$ appearing in the energy can now be replaced by the variables
$v^k$, leaving the only remaining dependence on the disordered $V$ in the
contribution of \eqref{eq:Vv.delta}, which is linear. The average over the
disorder can then be computed, which yields
\begin{equation}
  \begin{aligned}
    &\overline{\sum_{k=1}^M\sum_{a=1}^n\exp\left[i\int d1\,\hat v^k_a(1)V^k(\pmb\phi_a(1))\right]}
    \\
    &
    =\exp\left[
      -\frac12\sum_{k=1}^M\sum_{ab=1}^n\int d1\,d2\,\hat v_a^k(1)f\left(\frac{\pmb\phi_a(1)^T\pmb\phi_b(2)}N\right)\hat v_b^k(2)
    \right]
  \end{aligned}
\end{equation}
The result is factorized in the indices $k$ and Gaussian in the superfields $v$
and $\hat v$ with kernel
\begin{equation}
  \begin{bmatrix}
    B_a(1)\delta_{ab}\delta(1,2) & i\delta_{ab}\delta(1,2) \\
    i\delta_{ab}\delta(1,2) & f\left(\frac{\pmb\phi_a(1)^T\pmb\phi_b(2)}N\right)
  \end{bmatrix}
\end{equation}
Making the $M$ independent Gaussian integrals, we therefore have
\begin{equation}
  \begin{aligned}
    &\mathcal N(E,\mu)^n
    =\int\left(\prod_{a=1}^nd\hat\beta_a\,d\pmb\phi_a\right)
    \exp\bigg[
      \sum_a^n\hat\beta_aE \\
    &\qquad-\frac M2\log\operatorname{sdet}\left(
        \delta_{ab}\delta(1,2)+B_a(1)f\left(\frac{\pmb\phi_a(1)^T\pmb\phi_b(2)}N\right)
      \right)
    \bigg]
  \end{aligned}
\end{equation}
We make a change of variables from the fields $\pmb\phi$ to matrices $\mathbb Q_{ab}(1,2)=\frac1N\pmb\phi_a(1)^T\pmb\phi_b(2)$. This transformation results in a change of measure of the form
\begin{equation}
  \prod_{a=1}^n d\pmb\phi_a=d\mathbb Q\,(\operatorname{sdet}\mathbb Q)^\frac N2
  =d\mathbb Q\,\exp\left[\frac N2\log\operatorname{sdet}\mathbb Q\right]
\end{equation}
We therefore have
\begin{equation}
  \begin{aligned}
    &\mathcal N(E,\mu)^n
    =\int\left(\prod_{a=1}^nd\hat\beta_a\right)\,d\mathbb Q\,
    \exp\bigg[
      \sum_a^n\hat\beta_aE
      +\frac N2\log\operatorname{sdet}\mathbb Q
      \\
    &\qquad-\frac M2\log\operatorname{sdet}\left(
      \delta_{ab}\delta(1,2)+B_a(1)f(\mathbb Q_{ab}(1,2))
      \right)
    \bigg]
  \end{aligned}
\end{equation}
We now need to blow up our supermatrices into our physical order parameters. We have that
\begin{equation}
  \begin{aligned}
    &\mathbb Q_{ab}(1,2)
    =C_{ab}-G_{ab}(\bar\theta_1\theta_2+\bar\theta_2\theta_1) \\
    &\qquad-R_{ab}(\bar\theta_1\theta_1+\bar\theta_2\theta_2)
    -D_{ab}\bar\theta_1\theta_2\bar\theta_2\theta_2
  \end{aligned}
\end{equation}
where $C$, $R$, $D$, and $G$ are the matrices defined in
\eqref{eq:order.parameters}. Other possible combinations involving scalar
products between fermionic and bosonic variables do not contribute at physical
saddle points \cite{Kurchan_1992_Supersymmetry}. Inserting this expansion into
the expression above and evaluating the superdeterminants, we find
\begin{equation}
  \mathcal N(E,\mu)^n=\int d\hat\beta\,dC\,dR\,dD\,dG\,e^{nN\mathcal S_\mathrm{KR}(\hat\beta,C,R,D,G)}
\end{equation}
where the effective action is given by
\begin{widetext}
\begin{equation}
  \begin{aligned}
    &\mathcal S_\mathrm{KR}(\hat\beta,C,R,D,G)
    =\hat\beta E-\frac1n\operatorname{Tr}(G+R)\mu
    +\frac1n\frac12\Big(\log\det(CD+R^2)-\log\det G^2\Big)
    \\
    &-\frac1n\frac\alpha2\left\{\log\det\left[
      \Big(
        f'(C)\odot D-\hat\beta I+(G\odot G-R\odot R)\odot f''(C)
      \Big)f(C)
      +(I-R\odot f'(C))^2
    \right]-\log\det(I+G\odot f'(C))^2\right\}
  \end{aligned}
\end{equation}
where $\odot$ gives the Hadamard or componentwise product between the matrices, while other products and powers are matrix products and powers.

In the case where $\mu$ is not specified, the model has a BRST symmetry whose
Ward identities give $D=\hat\beta R$ and $G=-R$
\cite{Annibale_2004_Coexistence, Kent-Dobias_2023_How}. Using these relations,
the effective action becomes particularly simple:
\begin{equation}
  \mathcal S(\hat\beta, C, R)
  =
  \hat\beta E
  +\lim_{n\to0}\frac1n\left[
    -\frac\alpha2\log\det\left[
      I-\hat\beta f(C)(I-R\odot f'(C))^{-1}
    \right]
    +\frac12\log\det(I+\hat\beta CR^{-1})
  \right]
\end{equation}
This effective action is general for arbitrary matrices $C$ and $R$. When using
a replica symmetric ansatz of $C_{ab}=\delta_{ab}+c_0(1-\delta_{ab})$ and
$R_{ab}=r\delta_{ab}+r_0(1-\delta_{ab})$, the resulting function of
$\hat\beta$, $c_0$, $r$, and $r_0$ is
\begin{equation}
  \begin{aligned}
    \mathcal S=
    \hat\beta E
    -\frac\alpha 2\left[
      \log\left(1-\frac{\hat\beta\big(f(1)-f(c_0)\big)}{1-rf'(1)+r_0f'(c_0)}\right)
      -\frac{\hat\beta f(c_0)+r_0f'(c_0)}{
        1-\hat\beta\big(f(1)-f(c_0)\big)-rf'(1)+rf'(c_0)
      }+\frac{r_0f'(c_0)}{1-rf'(1)+r_0f'(c_0)}
    \right] \\
    +\frac12\left[
      \log\left(1+\frac{\hat\beta(1-c_0)}{r-r_0}\right)
      +\frac{\hat\beta c_0+r_0}{\hat\beta(1-c_0)+r-r_0}
      -\frac{r_0}{r-r_0}
    \right]
  \end{aligned}
\end{equation}
When $f(0)=0$ as in the cases directly studied in this work, this further
simplifies as $c_0=r_0=0$.

\end{widetext}

\bibliography{marginal}

\end{document}