#  Conditioning the complexity of random landscapes on marginal optima

Jaron Kent-Dobias<EMAIL_ADDRESS>Istituto Nazionale di
Fisica Nucleare, Sezione di Roma I, Rome, Italy 00184

###### 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 multispherical spin glasses, which are
Gaussian but non-GOE; and in sums of squared spherical random functions, which
are 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.

##  I Introduction

Systems with rugged landscapes are important across many disciplines, from the
physics of glasses and spin-glasses to statistical inference problems [1]. 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 [2, 3].

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. The level set associated with this threshold energy contains mostly
_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 [4, 5, 6]. 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 [7]. 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 be
used to bound out-of-equilibrium dynamical behavior.

Despite their importance in a wide variety of in and out of equilibrium
settings [8, 9, 10, 11, 12, 13, 14, 15, 16, 17], 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 [18]. However, this results in only a
characterization of the threshold energy and cannot characterize marginal
minima at other energies where they are a minority.

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 [19, 20]. 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 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 algorithms [21].

An outline of this paper follows. In Section II 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 GOE matrix with a shifted diagonal. In Section III we apply
this technique to the problem of characterizing marginal minima in random
landscapes. The following Section IV gives several examples of the marginal
complexity applied to specific models of increasing difficulty. Finally,
Section V summarizes this work and suggests necessary extensions.

##  II Conditioning on the smallest 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 II.2 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.

###  II.1 The general method

Consider an $N\times N$ real matrix $A$. An arbitrary function $g$ of the
minimum eigenvalue of $A$ can be expressed using integrals over
$\mathbf{s}\in\mathbb{R}^{N}$ as

$g(\lambda_{\textrm{min}}(A))=\lim_{\beta\to\infty}\int\frac{d\mathbf{s}\,\delta(N-\|\mathbf{s}\|^{2})e^{-\beta\mathbf{s}^{T}A\mathbf{s}}}{\int
d\mathbf{s}^{\prime}\,\delta(N-\|\mathbf{s}^{\prime}\|^{2})e^{-\beta\mathbf{s}^{\prime
T}A\mathbf{s}^{\prime}}}g\left(\frac{\mathbf{s}^{T}A\mathbf{s}}{N}\right)$ (1)

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

$\displaystyle\lim_{\beta\to\infty}\int\frac{d\mathbf{s}\,\delta(N-\|\mathbf{s}\|^{2})e^{-\beta\mathbf{s}^{T}A\mathbf{s}}}{\int
d\mathbf{s}^{\prime}\,\delta(N-\|\mathbf{s}^{\prime}\|^{2})e^{-\beta\mathbf{s}^{\prime
T}A\mathbf{s}^{\prime}}}g\left(\frac{\mathbf{s}^{T}A\mathbf{s}}{N}\right)$ (2)
$\displaystyle=\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}^{\prime}\,\delta(N-\|\mathbf{s}^{\prime}\|^{2})\mathbb{1}_{\operatorname{ker}(A-\lambda_{\mathrm{min}}(A)I)}(\mathbf{s}^{\prime})}g\left(\frac{\mathbf{s}^{T}A\mathbf{s}}{N}\right)$
$\displaystyle=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}^{\prime}\,\delta(N-\|\mathbf{s}^{\prime}\|^{2})\mathbb{1}_{\operatorname{ker}(A-\lambda_{\mathrm{min}}(A)I)}(\mathbf{s}^{\prime})}$
$\displaystyle=g(\lambda_{\mathrm{min}}(A))$

as desired. The first relation extends a technique for calculating the typical
minimum eigenvalue of an ensemble of matrices first introduced by Ikeda and
later used by Kent-Dobias in the context of random landscapes, and is similar
to an earlier technique for conditioning the value of the ground state energy
in random landscapes by Fyodorov and Le Doussal [22, 21, 24, 25]. 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}\|^{2}=N$ and the minimum
eigenspace,
$\mathbf{s}^{T}A\mathbf{s}=\mathbf{s}^{T}\mathbf{s}\lambda_{\mathrm{min}}(A)=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.

###  II.2 Simple example: shifted GOE

We demonstrate the efficacy of the technique by rederiving a well-known
result: the large-deviation function for pulling an eigenvalue from the bulk
of the GOE spectrum. Consider an ensemble of $N\times N$ matrices $A=B+\mu I$
for $B$ drawn from the GOE ensemble with entries whose variance is
$\sigma^{2}/N$. We know that the bulk spectrum of $A$ is a Wigner semicircle
with radius $2\sigma$ shifted by a constant $\mu$. Therefore, for
$\mu=2\sigma$, the minimum eigenvalue will typically be zero, while for
$\mu>2\sigma$ the minimum eigenvalue would need to be a large deviation from
the typical spectrum and its likelihood will be exponentially suppressed with
$N$. For $\mu<2\sigma$, the bulk of the typical spectrum contains zero and
therefore a larger $N^{2}$ deviation, moving an extensive number of
eigenvalues, would be necessary [26]. 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 large deviation function.

To compute this large deviation function, we will employ the method outlined
in the previous subsection to calculate

$\displaystyle e^{NG_{\lambda^{*}}(\mu)}$
$\displaystyle=P\big{(}\lambda_{\mathrm{min}}(B+\mu I)=\lambda^{*}\big{)}$ (3)
$\displaystyle=\overline{\delta\big{(}N\lambda^{*}-N\lambda_{\mathrm{min}}(B+\mu
I)\big{)}}$

where the overline is the average over $B$, and we have defined the large
deviation function $G_{\lambda^{*}}(\mu)$. Using the representation of
$\lambda_{\mathrm{min}}$ defined in (1), we have

$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}^{\prime}\,\delta(N-\|\mathbf{s}^{\prime}\|^{2})e^{-\beta\mathbf{s}^{\prime
T}(B+\mu
I)\mathbf{s}^{\prime}}}\,\delta\big{(}N\lambda^{*}-\mathbf{s}^{T}(B+\mu
I)\mathbf{s}\big{)}}$ (4)

Using replicas to treat the denominator ($x^{-1}=\lim_{m\to 0}x^{m-1}$) and
transforming the $\delta$-function to its Fourier representation, we have

$e^{NG_{\lambda^{*}}(\mu)}=\overline{\lim_{\beta\to\infty}\lim_{m\to 0}\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\\}}$ (5)

having introduced the auxiliary parameter $\hat{\lambda}$ in the Fourier
representation of the $\delta$-function. The whole expression, so transformed,
is an exponential integral linear in the matrix $B$. Taking the average over
$B$, we find

$\displaystyle e^{NG_{\lambda^{*}}(\mu)}=\lim_{\beta\to\infty}\lim_{m\to
0}\int
d\hat{\lambda}\prod_{\alpha=1}^{m}\left[d\mathbf{s}^{\alpha}\,\delta(N-\|\mathbf{s}^{\alpha}\|^{2})\right]$
(6) $\displaystyle\hskip
100.00015pt\times\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}^{2}N^{2}\right]\right\\}$

We make the Hubbard–Stratonovich transformation to the matrix field
$Q^{\alpha\beta}=\frac{1}{N}\mathbf{s}^{\alpha}\cdot\mathbf{s}^{\beta}$. This
produces an integral expression of the form

$e^{NG_{\lambda^{*}}(\mu)}=\lim_{\beta\to\infty}\lim_{m\to 0}\int
d\hat{\lambda}\,dQ\,e^{N\mathcal{U}_{\mathrm{GOE}}(\hat{\lambda},Q\mid\beta,\lambda^{*},\mu)}$
(7)

where the effective action $\mathcal{U}_{\mathrm{GOE}}$ is given by

$\displaystyle\mathcal{U}_{\textrm{GOE}}(\hat{\lambda},Q\mid\beta,\lambda^{*},\mu)=\hat{\lambda}(\lambda^{*}-\mu)+\lim_{m\to
0}\bigg{\\{}-m\beta\mu$ (8)
$\displaystyle+\sigma^{2}\left[\beta^{2}\sum_{\alpha\gamma}^{m}(Q^{\alpha\gamma})^{2}+2\beta\hat{\lambda}\sum_{\alpha}^{m}(Q^{1\alpha})^{2}+\hat{\lambda}^{2}\right]+\frac{1}{2}\log\det
Q\bigg{\\}}$

and $Q^{\alpha\alpha}=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 spherical model, but with the first
row singled out because of its unique coupling with $\hat{\lambda}$. The
resulting matrix has the form

$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}$ (9)

The relevant expressions in the effective action produce

$\displaystyle\sum_{\alpha\beta}(Q^{\alpha\beta})^{2}=m+2(m-1)\tilde{q}_{0}^{2}+(m-1)(m-2)q_{0}^{2}$
(10) $\displaystyle\sum_{\alpha}(Q^{1\alpha})^{2}=1+(m-1)\tilde{q}_{0}^{2}$
(11) $\displaystyle\begin{aligned} &\log\det Q=(m-2)\log(1-q_{0})\\\ &\hskip
60.00009pt+\log\big{[}1+(m-2)q_{0}-(m-1)\tilde{q}_{0}^{2}\big{]}\end{aligned}$
(12)

Inserting these expressions into the effective action and taking the limit of
$m$ to zero, we arrive at

$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\beta,\lambda^{*},\mu)}$
(13)

with the new effective action

$\displaystyle\mathcal{U}_{\mathrm{GOE}}(\hat{\lambda},q_{0},\tilde{q}_{0}\mid\beta,\lambda^{*},\mu)$
(14)
$\displaystyle\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]$
$\displaystyle\qquad-\log(1-q_{0})+\frac{1}{2}\log(1-2q_{0}+\tilde{q}_{0}^{2})$

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

$\displaystyle q_{0}$ $\displaystyle=1-y\beta^{-1}-z\beta^{-2}+O(\beta^{-3})$
(15) $\displaystyle\tilde{q}_{0}$
$\displaystyle=1-\tilde{y}\beta^{-1}-(z+\Delta z)\beta^{-2}+O(\beta^{-3})$
(16)

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

$\displaystyle\mathcal{U}_{\textrm{GOE}}(\hat{\lambda},y,\Delta
z\mid\infty,\lambda^{*},\mu)=\hat{\lambda}(\lambda^{*}-\mu)$ (17)
$\displaystyle\qquad+\sigma^{2}\big{[}\hat{\lambda}^{2}+4(y+\Delta
z)\big{]}+\frac{1}{2}\log\left(1-\frac{2\Delta z}{y^{2}}\right)$

Extremizing this action over the new parameters $y$, $\Delta z$, and
$\hat{\lambda}$, we find

$\displaystyle\hat{\lambda}$
$\displaystyle=\frac{1}{\sigma}\sqrt{\left(\frac{\mu-\lambda^{*}}{2\sigma}\right)^{2}-1}$
(18) $\displaystyle y$
$\displaystyle=\frac{1}{2\sigma}\left[\frac{\mu-\lambda^{*}}{2\sigma}+\sqrt{\left(\frac{\mu-\lambda^{*}}{2\sigma}\right)^{2}-1}\right]^{-1}$
(19) $\displaystyle\Delta z$
$\displaystyle=\frac{1}{4\sigma^{2}}\left[\left(\frac{\mu-\lambda^{*}}{2\sigma}\right)^{2}-1-\frac{\mu-\lambda^{*}}{2\sigma}\sqrt{\left(\frac{\mu-\lambda^{*}}{2\sigma}\right)^{2}-1}\right]$
(20)

Inserting this solution into the effective action we arrive at

$\displaystyle
G_{\lambda^{*}}(\mu)=\mathop{\textrm{extremum}}_{\hat{\lambda},y,\Delta
z}\mathcal{U}_{\mathrm{GOE}}(\hat{\lambda},y,\Delta
z\mid\infty,\lambda^{*},\mu)$ (21)
$\displaystyle=-\frac{\mu-\lambda^{*}}{2\sigma}\sqrt{\left(\frac{\mu-\lambda^{*}}{2\sigma}\right)^{2}-1}$
$\displaystyle\hskip
50.00008pt-\log\left[\frac{\mu-\lambda^{*}}{2\sigma}-\sqrt{\left(\frac{\mu-\lambda^{*}}{2\sigma}\right)^{2}-1}\right]$

This function is plotted in Fig. 1 for $\lambda^{*}=0$. For $\mu<2\sigma$,
$G_{0}(\mu)$ has an imaginary part. This indicates that the existence of a
zero minimum eigenvalue when $\mu<2\sigma$ 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 [26]. For $\mu\geq 2\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 typical bulk spectrum with zero, i.e.,
a pseudogap.

Figure 1:  The large deviation function $G_{0}(\mu)$ defined in (3) 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.

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 function will lie on the singular boundary between a purely
real and complex value.

###  II.3 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 marginality. 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=\frac{1}{N}\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, for a given $\mu$ the bottom of the spectrum can be
located by looking for the value $\lambda^{*}$ that maximizes the (real) large
deviation function. Inverting this reasoning, we can find the value
$\mu=\mu_{\textrm{m}}$ corresponding to a marginal spectrum by requiring that
the large deviation function has a maximum in $\lambda^{*}$ at
$\lambda^{*}=0$, or

$0=\frac{\partial}{\partial\lambda^{*}}G_{\lambda^{*}}(\mu_{\mathrm{m}})\bigg{|}_{\lambda^{*}=0}$
(22)

In the example problem of section II.2, 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$.

##  III Marginal complexity in random landscapes

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 III.1 we define the marginal complexity using the tools
of the previous section. In Section III.2 we review several general features
in a physicists’ approach to computing the marginal complexity. In Section
III.3 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.

###  III.1 Marginal complexity from Kac–Rice

The situation in the study of random landscapes is often as follows: an
ensemble of smooth energy functions $H:\mathbb{R}^{N}\to\mathbb{R}$ defines a
family of 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 typical geometry of landscapes drawn from
the ensemble 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 of identifying stationary points is reduced
to extremizing the Lagrangian

$L(\mathbf{x},\boldsymbol{\omega})=H(\mathbf{x})+\sum_{i=1}^{r}\omega_{i}g_{i}(\mathbf{x})$
(23)

with respect to $\mathbf{x}$ and the Lagrange multipliers
$\boldsymbol{\omega}=\\{\omega_{1},\ldots,\omega_{r}\\}$. The corresponding
gradient and Hessian of the energy associated with this constrained extremal
problem are

$\displaystyle\nabla H(\mathbf{x},\boldsymbol{\omega})=\partial
L(\mathbf{x},\boldsymbol{\omega})=\partial
H(\mathbf{x})+\sum_{i=1}^{r}\omega_{i}\partial g_{i}(\mathbf{x})$ (24)
$\displaystyle\begin{aligned}
\operatorname{Hess}H(\mathbf{x},\boldsymbol{\omega})&=\partial\partial
L(\mathbf{x},\boldsymbol{\omega})\\\ &=\partial\partial
H(\mathbf{x})+\sum_{i=1}^{r}\omega_{i}\partial\partial
g_{i}(\mathbf{x})\end{aligned}$ (25)

where $\partial=\frac{\partial}{\partial\mathbf{x}}$ will always represent the
derivative with respect to the vector argument $\mathbf{x}$.

The number of stationary points in a landscape for a particular function $H$
is found by integrating over the Kac–Rice measure

$\displaystyle d\nu_{H}(\mathbf{x},\boldsymbol{\omega})$ (26)
$\displaystyle\quad=d\mathbf{x}\,d\boldsymbol{\omega}\,\delta\big{(}\mathbf{g}(\mathbf{x})\big{)}\,\delta\big{(}\nabla
H(\mathbf{x},\boldsymbol{\omega})\big{)}\,\big{|}\det\operatorname{Hess}H(\mathbf{x},\boldsymbol{\omega})\big{|}$

with a $\delta$-function of the gradient and the constraints ensuring that we
count valid stationary points, and the determinant of the Hessian serving as
the Jacobian of the argument to the $\delta$ function [27, 28]. It is usually
more interesting to condition the count on interesting properties of the
stationary points, like the energy and spectrum trace, or

$\displaystyle d\nu_{H}(\mathbf{x},\boldsymbol{\omega}\mid E,\mu)$ (27)
$\displaystyle=d\nu_{H}(\mathbf{x},\boldsymbol{\omega})\,\delta\big{(}NE-H(\mathbf{x})\big{)}\,\delta\big{(}N\mu-\operatorname{Tr}\operatorname{Hess}H(\mathbf{x},\boldsymbol{\omega})\big{)}$

We specifically want to control the value of the minimum eigenvalue of the
Hessian at the stationary points. Using the method introduced in Section II,
we can write the number of stationary points with energy $E$, Hessian trace
$\mu$, and smallest eigenvalue $\lambda^{*}$ as

$\displaystyle\mathcal{N}_{H}(E,\mu,\lambda^{*})=\int
d\nu_{H}(\mathbf{x},\boldsymbol{\omega}\mid
E,\mu)\,\delta\big{(}N\lambda^{*}-\lambda_{\mathrm{min}}(\operatorname{Hess}H(\mathbf{x},\boldsymbol{\omega}))\big{)}$
(28) $\displaystyle=\lim_{\beta\to\infty}\int
d\nu_{H}(\mathbf{x},\boldsymbol{\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},\boldsymbol{\omega})\mathbf{s}}}{\int
d\mathbf{s}^{\prime}\,\delta(N-\|\mathbf{s}^{\prime}\|^{2})\delta(\mathbf{s}^{\prime
T}\partial\mathbf{g}(\mathbf{x}))e^{-\beta\mathbf{s}^{\prime
T}\operatorname{Hess}H(\mathbf{x},\boldsymbol{\omega})\mathbf{s}^{\prime}}}\delta\big{(}N\lambda^{*}-\mathbf{s}^{T}\operatorname{Hess}H(\mathbf{x},\boldsymbol{\omega})\mathbf{s}\big{)}$

where the additional $\delta$-functions

$\delta(\mathbf{s}^{T}\partial\mathbf{g}(\mathbf{x}))=\prod_{s=1}^{r}\delta(\mathbf{s}^{T}\partial
g_{i}(\mathbf{x}))$ (29)

ensure that the integrals involving potential eigenvectors $\mathbf{s}$ are
constrained to the tangent space of the configuration manifold at the point
$\mathbf{x}$.

The complexity of points with a specific energy, stability, and minimum
eigenvalue is defined as the average over the ensemble of functions $H$ of the
logarithm of the number $\mathcal{N}_{H}$ of stationary points, or

$\Sigma_{\lambda^{*}}(E,\mu)=\frac{1}{N}\overline{\log\mathcal{N}_{H}(E,\mu,\lambda^{*})}$
(30)

In practice, this can be computed by introducing replicas to treat the
logarithm ($\log x=\lim_{n\to 0}\frac{\partial}{\partial n}x^{n}$) and
introducing another set of replicas to treat each of the normalizations in the
numerator ($x^{-1}=\lim_{m\to-1}x^{m}$). This leads to the expression

$\displaystyle\Sigma_{\lambda^{*}}(E,\mu)$
$\displaystyle=\lim_{\beta\to\infty}\lim_{n\to
0}\frac{1}{N}\frac{\partial}{\partial
n}\int\prod_{a=1}^{n}\overline{\Bigg{[}d\nu_{H}(\mathbf{x}_{a},\boldsymbol{\omega}_{a}\mid
E,\mu)\,\delta\big{(}N\lambda^{*}-(\mathbf{s}_{a}^{1})^{T}\operatorname{Hess}H(\mathbf{x}_{a},\boldsymbol{\omega}_{a})\mathbf{s}_{a}^{1}\big{)}}$
(31) $\displaystyle\hskip 120.00018pt\overline{\times\lim_{m_{a}\to
0}\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},\boldsymbol{\omega}_{a})\mathbf{s}_{a}^{\alpha}}\right)\Bigg{]}}$

for the complexity of stationary points of a given energy, trace, and smallest
eigenvalue.

The marginal complexity follows from the complexity as a function of $\mu$ and
$\lambda^{*}$ in an analogous way to Section II.3. In general, one sets
$\lambda^{*}=0$ and tunes $\mu$ from a sufficiently large value until the
complexity develops an imaginary component, which corresponds to the bulk of
the spectrum touching zero. The value $\mu=\mu_{\mathrm{m}}$ that satisfies
this is the marginal stability.

In the cases studied here with zero signal-to-noise, a simpler approach is
possible. The marginal stability $\mu=\mu_{\text{m}}$ can be identified by
requiring that the complexity is stationary with respect to changes in the
value of the minimum eigenvalue $\lambda^{*}$, or

$0=\frac{\partial}{\partial\lambda^{*}}\Sigma_{\lambda^{*}}(E,\mu_{\text{m}}(E))\bigg{|}_{\lambda^{*}=0}$
(32)

The marginal complexity follows by evaluating the complexity conditioned on
$\lambda^{*}=0$ at the marginal stability $\mu=\mu_{\text{m}}(E)$,

$\Sigma_{\text{m}}(E)=\Sigma_{0}(E,\mu_{\text{m}}(E))$ (33)

###  III.2 General features of saddle point computation

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 components of the Kac–Rice measure defined in (26) and (27) into elements
of an exponential integral over configuration space. To begin with, all Dirac
$\delta$ functions are expressed using their Fourier representation, with

$\displaystyle\delta\big{(}\nabla
H(\mathbf{x}_{a},\boldsymbol{\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},\boldsymbol{\omega}_{a})}$ (34)
$\displaystyle\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}))}$
(35) $\displaystyle\begin{aligned}
&\delta\big{(}N\lambda^{*}-(\mathbf{s}_{a}^{1})^{T}\operatorname{Hess}H(\mathbf{x}_{a},\boldsymbol{\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},\boldsymbol{\omega})\mathbf{s}_{a}^{1})}\end{aligned}$
(36)

To do this we have introduced auxiliary fields $\hat{\mathbf{x}}_{a}$,
$\hat{\beta}_{a}$, and $\hat{\lambda}_{a}$. Because the permutation symmetry
of replica vectors is preserved in rsb orders, the order parameters
$\hat{\beta}$ and $\hat{\lambda}$ will quickly lose their indices, since they
will ubiquitously be constant over the replica index at the eventual saddle
point solution.

We would like to make a similar treatment of the determinant of the Hessian
that appears in (26). The standard approach is to drop the absolute value
function around the determinant. This can potentially lead to severe problems
with the complexity [19]. However, it is a justified step when the parameters
of the problem $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 the present problem of marginal minima, which lie right at
the edge of disaster.

Dropping the absolute value function allows us to write

$\det\operatorname{Hess}H(\mathbf{x}_{a},\boldsymbol{\omega}_{a})=\int
d\bar{\boldsymbol{\eta}}_{a}\,d\boldsymbol{\eta}_{a}\,e^{-\bar{\boldsymbol{\eta}}_{a}^{T}\operatorname{Hess}H(\mathbf{x}_{a},\boldsymbol{\omega}_{a})\boldsymbol{\eta}_{a}}$
(37)

using $N$-dimensional Grassmann vectors $\bar{\boldsymbol{\eta}}_{a}$ and
$\boldsymbol{\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 [4, 7]. 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
the trace of $\partial\partial g_{i}$. The result is that

$\mu=\frac{1}{N}\operatorname{Tr}\operatorname{Hess}H(\mathbf{x})=\frac{1}{N}\sum_{i=1}^{r}\omega_{i}\operatorname{Tr}\partial\partial
g_{i}(\mathbf{x})+O(N^{-1})$ (38)

In particular, here we study only cases with quadratic $g_{i}$, which results
in a linear 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.

###  III.3 Superspace representation

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 [29, 30, 31]. Appendix A introduces the notation and methods of
superspace algebra. Here we describe how it can be used to simplify the
complexity calculation for marginal minima.

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

$\boldsymbol{\phi}_{a}^{\alpha}(1,2)=\mathbf{x}_{a}+\bar{\theta}_{1}\boldsymbol{\eta}_{a}+\bar{\boldsymbol{\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})$
(39)

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

$\displaystyle\mathcal{N}_{H}(E,\mu,\lambda^{*})^{n}$
$\displaystyle=\lim_{\beta\to\infty}\int
d\boldsymbol{\omega}\,d\hat{\beta}\,d\hat{\lambda}\prod_{a=1}^{n}\lim_{m_{a}\to
0}\prod_{\alpha=1}^{m_{a}}d\boldsymbol{\phi}_{a}^{\alpha}\exp\left\\{\delta^{\alpha
1}N(\hat{\beta}E+\hat{\lambda}\lambda^{*})+\int
d1\,d2\,B^{\alpha}(1,2)L\big{(}\boldsymbol{\phi}_{a}^{\alpha}(1,2),\boldsymbol{\omega}\big{)}\right\\}$
(40)

Here we have also defined the operator

$B^{\alpha}(1,2)=\delta^{\alpha
1}\bar{\theta}_{2}\theta_{2}(1-\hat{\beta}\bar{\theta}_{1}\theta_{1})-\delta^{\alpha
1}\hat{\lambda}-\beta$ (41)

which encodes various aspects of the complexity problem. When the Lagrangian
is expanded in a series with respect to the Grassmann indices and the
definition of $B$ inserted, the result of the Grassmann integrals produces
exactly the content of the integrand in (31) with the substitutions (34),
(35), (36), and (37) of the Dirac $\delta$ functions and the determinant made.
The new measures

$\displaystyle d\boldsymbol{\phi}_{a}^{\alpha}$
$\displaystyle=\left[d\mathbf{x}_{a}\,\delta\big{(}\mathbf{g}(\mathbf{x}_{a})\big{)}\,\frac{d\hat{\mathbf{x}}_{a}}{(2\pi)^{N}}\,d\boldsymbol{\eta}_{a}\,d\bar{\boldsymbol{\eta}}_{a}\,\delta^{\alpha
1}+(1-\delta^{\alpha
1})\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{)}$
(42) $\displaystyle d\boldsymbol{\omega}$
$\displaystyle=\bigg{(}\prod_{i=1}^{r}d\omega_{i}\bigg{)}\,\delta\bigg{(}N\mu-\sum_{i}^{r}\omega_{i}\operatorname{Tr}\partial\partial
g_{i}\bigg{)}$ (43)

collect the individual measures of the various fields embedded in the
superfield, along with their constraints. With this way of writing the
replicated count, the problem of marginal complexity temporarily takes the
schematic form of an equilibrium calculation with configurations
$\boldsymbol{\phi}$, inverse temperature $B$, and energy $L$. This makes the
intermediate pieces of the calculation dramatically simpler. Of course the
intricacies of the underlying problem are not banished: near the end of the
calculation, terms involving the superspace must be expanded. We will make use
of this representation to simplify the analysis of the marginal complexity
when analyzing random sums of squares in Section IV.3.

##  IV Examples

In this section we present analysis of marginal complexity in three random
landscapes. In Section IV.1 we treat 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 IV.2 we apply the methods to a
multispherical spin glass, which is still Gaussian but has a non-GOE spectrum
whose shape can vary between stationary points. Finally, in Section IV.3 we
analyze a model of sums of squared random functions, which is non-Gaussian and
whose Hessian statistics depend on the conditioning of the energy and
gradient.

###  IV.1 Spherical spin glasses

The spherical spin glasses are a family of models that encompass every
isotropic Gaussian field on the hypersphere. Their configuration space is the
sphere $S^{N-1}$ defined by all $\mathbf{x}\in\mathbb{R}^{N}$ such that
$0=g(\mathbf{x})=\frac{1}{2}(\|\mathbf{x}\|^{2}-N)$. One can consider the
models as defined by ensembles of centered Gaussian functions $H$ such that
the covariance between two points in the configuration space is

$\overline{H(\mathbf{x})H(\mathbf{x}^{\prime})}=Nf\left(\frac{\mathbf{x}\cdot\mathbf{x}^{\prime}}{N}\right)$
(44)

for some function $f$ with positive series coefficients. Such functions can be
considered to be made up of all-to-all tensorial interactions, with

$H(\mathbf{x})=\sum_{p=0}^{\infty}\frac{1}{p!}\sqrt{\frac{f^{(p)}(0)}{N^{p-1}}}\sum_{i_{1}\cdots
i_{p}}^{N}J_{i_{1}\cdots i_{p}}x_{i_{1}}\cdots x_{i_{p}}$ (45)

and the elements of the tensors $J$ being independently distributed with the
unit normal distribution [32]. We focus on marginal minima in models with
$f^{\prime}(0)=0$, which corresponds to models without a random external
field. Such a random field would correspond in each individual sample $H$ to a
signal, and therefore complicate the analysis by correlating the positions of
stationary points and the eigenvectors of their Hessians.

The marginal optima of these models can be studied without the methods
introduced in this paper, and have been in the past [4, 7]. First, these
models are Gaussian, so at large $N$ the Hessian is statistically independent
of the gradient and energy [19, 20]. 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^{\prime\prime}(1)}$. Therefore, all marginal minima
in these systems have the same constant shift $\mu=\mu_{\mathrm{m}}$. Despite
the fact that 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 some things about its application in more nontrivial
settings.

Note that in the pure version of these models with $f(q)=\frac{1}{2}q^{p}$,
the methods of this section must be amended slightly. This is because in these
models there is an exact correspondence $\mu=-pE$ between the trace of the
Hessian and the energy, and therefore they cannot be fixed independently. This
correspondence implies that when $\mu=\mu_{\mathrm{m}}$, the corresponding
energy level $E_{\mathrm{th}}=-\frac{1}{p}\mu_{\mathrm{m}}$ contains all
marginal minima. This is what gives this threshold energy such singular
importance to dynamics in the pure spherical models.

The procedure to treat the complexity of the spherical models has been made in
detail elsewhere [7]. Here we make only a sketch of the steps involved. First
we notice that $\mu=\frac{1}{N}\omega\operatorname{Tr}\partial\partial
g(\mathbf{x})=\omega$, so that the only Lagrange multiplier $\omega$ in this
problem is set directly to the shift $\mu$. The substitutions (34), (35), and
(36) are made to convert the Dirac $\delta$ functions into exponential
integrals, and the substitution (37) is made to likewise convert the
determinant.

Once these substitutions have been made, the entire expression (31) 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 associated with replica $a$ into the linear functional
$\mathcal{O}_{a}$ then the average over the ensemble of functions $H$ gives

$\displaystyle\overline{e^{\sum_{a}^{n}\mathcal{O}_{a}H(\mathbf{x}_{a})}}$
$\displaystyle=e^{\frac{1}{2}\sum_{a}^{n}\sum_{b}^{n}\mathcal{O}_{a}\mathcal{O}_{b}\overline{H(\mathbf{x}_{a})H(\mathbf{x}_{b})}}$
(46)
$\displaystyle=e^{N\frac{1}{2}\sum_{a}^{n}\sum_{b}^{n}\mathcal{O}_{a}\mathcal{O}_{b}f\big{(}\frac{\mathbf{x}_{a}\cdot\mathbf{x}_{b}}{N}\big{)}}$

The result is an integrand that depends on the many vector variables we have
introduced only 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

$\displaystyle C_{ab}=\frac{1}{N}\mathbf{x}_{a}\cdot\mathbf{x}_{b}\quad$
$\displaystyle
R_{ab}=-i\frac{1}{N}\mathbf{x}_{a}\cdot\hat{\mathbf{x}}_{b}\quad$
$\displaystyle
D_{ab}=\frac{1}{N}\hat{\mathbf{x}}_{a}\cdot\hat{\mathbf{x}}_{b}$
$\displaystyle
Q_{ab}^{\alpha\gamma}=\frac{1}{N}\mathbf{s}_{a}^{\alpha}\cdot\mathbf{s}_{b}^{\gamma}\quad$
$\displaystyle\hat{X}^{\alpha}_{ab}=-i\frac{1}{N}\hat{\mathbf{x}}_{a}\cdot\mathbf{s}_{b}^{\alpha}\quad$
$\displaystyle
X^{\alpha}_{ab}=\frac{1}{N}\mathbf{x}_{a}\cdot\mathbf{s}_{b}^{\alpha}$
$\displaystyle
G_{ab}=\frac{1}{N}\bar{\boldsymbol{\eta}}_{a}\cdot\boldsymbol{\eta}_{b}$ (47)

Order parameters that mix the normal and Grassmann variables generically
vanish in these settings and we don’t consider them here [33]. This
transformation changes the measure of the integral, with

$\displaystyle\prod_{a=1}^{n}d\mathbf{x}_{a}\,\frac{d\hat{\mathbf{x}}_{a}}{(2\pi)^{N}}\,d\bar{\boldsymbol{\eta}}_{a}\,d\boldsymbol{\eta}_{a}\,\prod_{\alpha=1}^{m_{a}}d\mathbf{s}_{a}^{\alpha}$
(48) $\displaystyle\quad=dC\,dR\,dD\,dG\,dQ\,dX\,d\hat{X}\,(\det J)^{N/2}(\det
G)^{-N}$

where $J$ is the Jacobian of the transformation in the real-valued fields.
This Jacobian takes a block form

$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}$ (49)

The Grassmann integrals produces their own inverted Jacobian. The matrix that
make up the blocks of the matrix $J$ are such that $C$, $R$, and $D$ are
$n\times n$ matrices indexed by their lower indices, $Q_{ab}$ is an
$m_{a}\times m_{b}$ matrix indexed by its upper indices, while $X_{a}$ is an
$n\times m_{a}$ matrix with one lower and one upper index.

These steps follow identically to those more carefully outlined in the cited
papers [4, 7]. Following them in the present case, we arrive at a form for the
complexity of stationary points with fixed energy $E$, stability $\mu$, and
lowest eigenvalue $\lambda^{*}$ with

$\displaystyle\Sigma_{\lambda^{*}}(E,\mu)=\lim_{\beta\to\infty}\lim_{n\to
0}\lim_{m_{1}\cdots m_{n}\to 0}\frac{1}{N}\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)$ (50)
$\displaystyle\qquad+nN\mathcal{U}_{\mathrm{SSG}}(\hat{\lambda},Q,X,\hat{X}\mid\beta,\lambda^{*},\mu,C)+\frac{N}{2}\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{\\}}$

The exponential integrand is split into two effective actions coupled only by
a residual determinant. The first of these actions is the usual effective
action for the complexity of the spherical spin glasses, or

$\displaystyle\mathcal{S}_{\mathrm{SSG}}(\hat{\beta},C,R,D,G\mid
E,\mu)=\hat{\beta}E+\lim_{n\to
0}\frac{1}{n}\bigg{\\{}-\mu\operatorname{Tr}(R+G)$ (51)
$\displaystyle\qquad+\frac{1}{2}\sum_{ab}\left(\hat{\beta}^{2}f(C_{ab})+\big{(}2\hat{\beta}R_{ab}-D_{ab}\big{)}f^{\prime}(C_{ab})+(R_{ab}^{2}-G_{ab}^{2})f^{\prime\prime}(C_{ab})\right)+\frac{1}{2}\log\det\begin{bmatrix}C&iR\\\
iR^{T}&D\end{bmatrix}-\log\det G\bigg{\\}}$

The second of these actions is analogous to the effective action (8) from the
GOE example of Section II.2 and contains the contributions from the marginal
pieces of the calculation, and is given by

$\displaystyle\mathcal{U}_{\mathrm{SSG}}(\hat{\lambda},Q,X,\hat{X}\mid\beta,\lambda^{*},\mu,C)=\hat{\lambda}\lambda^{*}+\lim_{n\to
0}\lim_{m_{1}\cdots m_{n}\to 0}\frac{1}{n}\Bigg{\\{}\frac{1}{2}\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}^{n}f^{\prime\prime}(C_{ab})$ (52)
$\displaystyle\quad\times\Bigg{[}\beta\sum_{\alpha}^{m_{a}}\left(\beta\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}^{\alpha
1}+\sum_{\alpha}^{m_{b}}Q_{ab}^{1\alpha}\right)\Bigg{]}\Bigg{\\}}$

The fact that the complexity can be split into two relatively independent
pieces in this way is a characteristic of the Gaussian nature of the spherical
spin glass. In Section IV.3 we will study a model whose energy is not Gaussian
and where such a decomposition is impossible.

There are some dramatic simplifications that emerge from the structure of this
particular problem. First, notice that 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 without a fixed or
random field, 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 a landscape without a signal, where no direction is any more important than
any other, we expect such correlations to be zero: where a state is located
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 _and_
their soft eigenvectors. Therefore, in these instances one must account for
solutions with nonzero $X$ and $\hat{X}$.

We similarly expect that $Q_{ab}=0$ for $a\neq b$. For the contrary to be
true, eigenvectors at independently sampled stationary points would need to
have their directions correlated. This is expected in situations with a
signal, where such correlations would be driven by a shared directional bias
towards the signal. In the present situation, where there is no signal, such
correlations do not exist.

When we take $X=\hat{X}=0$ and
$Q^{\alpha\beta}_{ab}=\delta_{ab}Q^{\alpha\beta}$, we find that

$\mathcal{U}_{\mathrm{SSG}}(\hat{\lambda},Q,0,0\mid\beta,\lambda^{*},\mu,C)=\mathcal{U}_{\mathrm{GOE}}(\hat{\lambda},Q\mid\beta,\lambda^{*},\mu)$
(53)

with $\sigma^{2}=f^{\prime\prime}(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 GOE problem. This is perhaps not so
surprising, since we established from the beginning that the Hessian of the
spherical spin glasses belongs to the GOE class.

The remaining analysis of the eigenvalue-dependent part
$\mathcal{U}_{\mathrm{SSG}}$ follows precisely the same steps as were made in
Section II.2 for the GOE example. The result of the calculation is also the
same: the exponential factor containing $\mathcal{U}_{\mathrm{SSG}}$ produces
precisely the large deviation function $G_{\lambda^{*}}(\mu)$ of (21) (again
with $\sigma^{2}=f^{\prime\prime}(1)$). The remainder of the integrand
depending on $\mathcal{S}_{\mathrm{SSG}}$ produces the ordinary complexity of
the spherical spin glasses without conditions on the Hessian eigenvalue. We
therefore find that

$\Sigma_{\lambda^{*}}(E,\mu)=\Sigma(E,\mu)+G_{\lambda^{*}}(\mu)$ (54)

We find the marginal complexity by solving

$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}$
(55)

which gives $\mu_{\mathrm{m}}(E)=2\sigma=2\sqrt{f^{\prime\prime}(1)}$
independent of $E$, as we presaged above. Since $G_{0}(\mu_{\mathrm{m}})=0$,
this gives finally

$\Sigma_{\mathrm{m}}(E)=\Sigma_{0}(E,\mu_{\mathrm{m}}(E))=\Sigma(E,\mu_{\mathrm{m}})$
(56)

that the marginal complexity in these models is simply the ordinary complexity
evaluated at a fixed trace $\mu_{\mathrm{m}}$ of the Hessian.

###  IV.2 Multispherical spin glasses

The multispherical spin glasses 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}$. The
two spheres give rise to two constraints: for
$\mathbf{x}=[\mathbf{x}^{(1)},\mathbf{x}^{(2)}]$ with components
$\mathbf{x}^{(1)},\mathbf{x}^{(2)}\in\mathbb{R}^{N}$, the constraints are
$0=g_{1}(\mathbf{x})=\frac{1}{2}(\|\mathbf{x}^{(1)}\|^{2}-N)$ and
$0=g_{2}(\mathbf{x})=\frac{1}{2}(\|\mathbf{x}^{(2)}\|^{2}-N)$. These two
constraints are fixed by two Lagrange multipliers $\omega_{1}$ and
$\omega_{2}$.

The energy in our multispherical spin glass is given by

$H(\mathbf{x})=H_{1}(\mathbf{x}^{(1)})+H_{2}(\mathbf{x}^{(2)})-\epsilon\mathbf{x}^{(1)}\cdot\mathbf{x}^{(2)}$
(57)

The energy $H_{i}$ of each individual sphere is taken to be a centered
Gaussian random function with a covariance given in the usual spherical spin
glass way for $\mathbf{x},\mathbf{x}^{\prime}\in\mathbb{R}^{N}$ by

$\overline{H_{i}(\mathbf{x})H_{j}(\mathbf{x}^{\prime})}=N\delta_{ij}f_{i}\left(\frac{\mathbf{x}\cdot\mathbf{x}^{\prime}}{N}\right)$
(58)

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 [34, 35,
36, 37, 38, 39, 40]. As far as we are aware, the deterministically coupled
model has not been previously studied, except as a thought experiment in [7].

We again make use of the method of Lagrange multipliers to find stationary
points on the constrained configuration space. The Lagrangian and its gradient
and Hessian are

$\displaystyle\begin{aligned}
L(\mathbf{x})&=H(\mathbf{x})+\frac{1}{2}\omega_{1}\big{(}\|\mathbf{x}^{(1)}\|^{2}-N\big{)}\\\
&\qquad\qquad\qquad+\frac{1}{2}\omega_{2}\big{(}\|\mathbf{x}^{(2)}\|^{2}-N\big{)}\end{aligned}$
(59) $\displaystyle\nabla
H(\mathbf{x},\boldsymbol{\omega})=\begin{bmatrix}\partial_{1}H_{1}(\mathbf{x}^{(1)})-\epsilon\mathbf{x}^{(2)}+\omega_{1}\mathbf{x}^{(1)}\\\
\partial_{2}H_{2}(\mathbf{x}^{(2)})-\epsilon\mathbf{x}^{(1)}+\omega_{2}\mathbf{x}^{(2)}\end{bmatrix}$
(60) $\displaystyle\begin{aligned}
&\operatorname{Hess}H(\mathbf{x},\boldsymbol{\omega})\\\
&\quad=\begin{bmatrix}\partial_{1}\partial_{1}H_{1}(\mathbf{x}^{(1)})+\omega_{1}I&-\epsilon
I\\\ -\epsilon
I&\partial_{2}\partial_{2}H_{2}(\mathbf{x}^{(2)})+\omega_{2}I\end{bmatrix}\end{aligned}$
(61)

where $\partial_{1}=\frac{\partial}{\partial\mathbf{x}^{(1)}}$ and
$\partial_{2}=\frac{\partial}{\partial\mathbf{x}^{(2)}}$. Like in the
spherical spin glasses, 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 fixing $\mu$.

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 C 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 C and find 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 single 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

$Q^{ij,\alpha\gamma}_{ab}=\frac{1}{N}\mathbf{s}^{(i),\alpha}_{a}\cdot\mathbf{s}^{(j),\gamma}_{b}$
(62)

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. We
immediately make 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 for the complexity has the schematic form of
(50), but with different effective actions depending now on overlaps inside
each of the two spheres and between the two spheres. The effective action for
the traditional complexity of the multispherical spin glass is

$\displaystyle\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},G^{21}\mid
E,\omega_{1},\omega_{2})=\hat{\beta}(E-E_{1}-E_{2}-\epsilon c_{d}^{12})$
$\displaystyle+\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})+\lim_{n\to
0}\frac{1}{n}\Bigg{\\{}\epsilon\operatorname{Tr}(R^{12}+R^{21}+G^{12}+G^{21}-\hat{\beta}C^{12})$
$\displaystyle\quad+\frac{1}{2}\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})\Bigg{\\}}$
(63)

which is the sum of two effective actions (51) for the spherical spin glass
associated with each individual sphere, and some coupling terms. The order
parameters are defined the same as in the spherical spin glasses, but now with
raised indices to indicate whether the vectors come from one or the other
spherical subspace. The effective action for the eigenvalue-dependent part of
the complexity is likewise given by

$\displaystyle\mathcal{U}_{\mathrm{MSG}}(\hat{q},\hat{\lambda},Q^{11},Q^{22},Q^{12}\mid\beta,\lambda^{*},\omega_{1},\omega_{2})$
$\displaystyle\quad=\lim_{m\to
0}\bigg{\\{}\sum_{\alpha=1}^{m}\left[\hat{q}^{\alpha}(Q^{11,\alpha\alpha}+Q^{22,\alpha\alpha}-1)-\beta(\omega_{1}Q^{11,\alpha\alpha}+\omega_{2}Q^{22,\alpha\alpha}-2\epsilon
Q^{12,\alpha\alpha})\right]-\hat{\lambda}(\omega_{1}Q^{11,11}+\omega_{2}Q^{22,11}-2\epsilon
Q^{12,11})$
$\displaystyle\qquad\qquad+\sum_{i=1,2}f_{i}^{\prime\prime}(1)\left[\beta^{2}\sum_{\alpha\gamma}^{m}(Q^{ii,\alpha\gamma})^{2}+2\beta\hat{\lambda}\sum_{\alpha}^{m}(Q^{ii,1\alpha})^{2}+\hat{\lambda}^{2}(Q^{ii,11})^{2}\right]+\frac{1}{2}\log\det\begin{bmatrix}Q^{11}&Q^{12}\\\
Q^{12}&Q^{22}\end{bmatrix}\bigg{\\}}$ (64)

The new variables $\hat{q}^{\alpha}$ are Lagrange multipliers introduced to
enforce the constraint that $Q^{11,\alpha\alpha}+Q^{22,\alpha\alpha}=1$.
Because of this constraint, the diagonal of the $Q$ matrices cannot be taken
to be 1 as in Section II.2. Instead we take each of the matrices $Q^{11}$,
$Q^{22}$, and $Q^{12}$ to have the planted replica symmetric form of (9), but
with the diagonal not necessarily equal to 1, so

$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}$ (65)

This requires us to introduce two new order parameters $\tilde{q}^{ij}_{d}$
and $q^{ij}_{d}$ per pair $(i,j)$, in addition to the off-diagonal order
parameters $\tilde{q}_{0}^{ij}$ and $q_{0}^{ij}$ already present in (9). We
also need two separate Lagrange multipliers $\hat{q}$ and $\hat{\tilde{q}}$ to
enforce the tangent space normalization $q_{d}^{11}+q_{d}^{22}=1$ and
$\tilde{q}_{d}^{11}+\tilde{q}_{d}^{22}=1$ for the tilde and untilde replicas,
respectively, which will in general take different values at the saddle point.
When this ansatz is inserted into the expression (64) for the effective action
and the limit of $m\to 0$ is taken, we find

$\displaystyle\mathcal{U}_{\mathrm{MSG}}(\hat{q},\hat{\tilde{q}},\hat{\lambda},\tilde{q}_{d}^{11},\tilde{q}_{0}^{11},q_{d}^{11},q_{0}^{11},\tilde{q}_{d}^{22},\tilde{q}_{0}^{22},q_{d}^{22},q_{0}^{22},\tilde{q}_{d}^{12},\tilde{q}_{0}^{12},q_{d}^{12},q_{0}^{12}\mid\beta,\lambda^{*},\omega_{1},\omega_{2})$
$\displaystyle=\sum_{i=1,2}\left\\{f_{i}^{\prime\prime}(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}\tilde{q}^{ii}_{d}\omega_{i}-\beta(\tilde{q}^{ii}_{d}-q^{ii}_{d})\omega_{i}\right\\}$
$\displaystyle+\frac{1}{2}\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}_{0}q^{22}_{0}+\tilde{q}^{11}_{d}\tilde{q}^{22}_{0}+\tilde{q}^{11}_{0}q^{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)$
$\displaystyle\qquad\qquad+2\left(3(q^{12}_{0})^{2}-(\tilde{q}^{12}_{0})^{2}-2q^{12}_{0}q^{12}_{d}-3q^{11}_{0}q^{22}_{0}+q^{11}_{d}q^{22}_{0}+\tilde{q}^{11}_{0}\tilde{q}^{22}_{0}+q^{11}_{0}q^{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)$
$\displaystyle\qquad\qquad-\left(2(q^{12}_{0})^{2}-(\tilde{q}^{12}_{0})^{2}-(q^{12}_{d})^{2}-2q^{11}_{0}q^{22}_{0}+\tilde{q}^{11}_{0}\tilde{q}^{22}_{0}+q^{11}_{d}q^{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{]}$
$\displaystyle-\log\left[(q^{11}_{d}-q^{11}_{0})(q^{22}_{d}-q^{22}_{0})-(q^{12}_{d}-q^{12}_{0})^{2}\right]+2\epsilon\big{[}\hat{\lambda}\tilde{q}^{12}_{d}+\beta(\tilde{q}^{12}_{d}-q^{12}_{d})\big{]}-\hat{q}(q^{11}_{d}+q^{22}_{d}-1)+\hat{\tilde{q}}(\tilde{q}^{11}_{d}+\tilde{q}^{22}_{d}-1)$
(66)

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}$. Notice that
in this case, the asymptotic behavior of the off-diagonal elements is to
approach the value of the diagonal rather than to approach 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 terms also approache the same diagonal value as
the untilde terms, but with potentially different rates.

As before, in order for the logarithmic term to stay finite, there are
necessary constraints on the values $y$. These are

$\displaystyle\frac{1}{2}(y^{11}_{d}-\tilde{y}^{11}_{d})=y^{11}_{0}-\tilde{y}^{11}_{0}$
(67)
$\displaystyle\frac{1}{2}(y^{22}_{d}-\tilde{y}^{22}_{d})=y^{22}_{0}-\tilde{y}^{22}_{0}$
(68)
$\displaystyle\frac{1}{2}(y^{12}_{d}-\tilde{y}^{12}_{d})=y^{12}_{0}-\tilde{y}^{12}_{0}$
(69)

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.

When the $\beta$-dependence of the $q$ variables is inserted into the
effective action (66) and the limit $\beta\to\infty$ taken, we find an
expression that is too large to report here. However, it can be extremized
over all of the variables in the problem just as in the previous examples to
find the values of the Lagrange multipliers $\omega_{1}$ and $\omega_{2}$
corresponding to marginal minima. Fig. 2(a) shows examples of the $\omega_{1}$
and $\omega_{2}$ corresponding to marginal spectra for a variety of couplings
$\epsilon$ when the covariances of the energy on the two spherical subspaces
are such that $1=f_{1}^{\prime\prime}(1)=f_{2}^{\prime\prime}(1)$. Fig. 2(b)
shows the Hessian spectra associated with some specific pairs
$(\omega_{1},\omega_{2})$. When $\epsilon=0$ and the two spheres are
uncoupled, we find the result for two independent spherical spin glasses: if
either $\omega_{1}=2\sqrt{f^{\prime\prime}(1)}=2$ or
$\omega_{2}=2\sqrt{f^{\prime\prime}(1)}=2$ and the other Lagrange multiplier
is larger than 2, then we have a marginal minimum made up of the Cartesian
product of a marginal minimum on one subspace and a stable minimum on the
other.

Figure 2:  Properties of marginal minima in the multispherical model. (a)
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}^{\prime\prime}(1)=1$,
$\sigma_{2}^{2}=f_{2}^{\prime\prime}(1)=1$, and various $\epsilon$. (b)
Spectra corresponding to the parameters $\omega_{1}$ and $\omega_{2}$ marked
by the circles in panel (a). (c) The complexity of marginal minima in a
multispherical model with $f_{1}(q)=\frac{1}{6}q^{3}$ and
$f_{2}(q)=\frac{1}{12}q^{4}$ for a variety of $\epsilon$. Since
$f_{1}^{\prime\prime}(1)=f_{2}^{\prime\prime}(1)=1$, the marginal values
correspond precisely to those in (a–b).

Fig. 2(c) shows the complexity of marginal minima in an example where both
$H_{1}$ and $H_{2}$ correspond to pure $p$-spin models, with
$f_{1}(q)=\frac{1}{6}q^{3}$ and $f_{2}(q)=\frac{1}{12}q^{4}$. Despite having
different covariance functions, these both satisfy
$1=f_{1}^{\prime\prime}(1)=f_{2}^{\prime\prime}(1)$ and therefore have
marginal minima for Lagrange multipliers that satisfy the relationships in
Fig.2(a). In the uncoupled system with $\epsilon=0$, the most common type of
marginal stationary point consists of independently marginal stationary points
in the two subsystems, with $\omega_{1}=\omega_{2}=2$. As $\epsilon$ is
increased, the most common type of marginal minimum drifts toward points with
$\omega_{1}>\omega_{2}$.

Multispherical spin glasses may be an interesting platform for testing ideas
about which among the possible marginal minima can dynamics, and cannot. In
the limit where $\epsilon=0$ and the configurations of the two spheres are
independent, the minima found dynamically should be marginal on both
subspaces. Just because technically on the expanded configuration space the
Cartesian product of a deep 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.

###  IV.3 Sums of squared random functions

In this subsection we consider perhaps the simplest example of a non-Gaussian
landscape: the problem of sums of squared random functions. This problem has a
close resemblance to nonlinear least squares optimization. Though, for reasons
we will see it is easier to make predictions for nonlinear _most_ squares,
i.e., the problem of maximizing the sum of squared terms. We again take a
spherical configuration space with $\mathbf{x}\in S^{N-1}$ and
$0=g(\mathbf{x})=\frac{1}{2}(\|\mathbf{x}\|^{2}-N)$ as in the spherical spin
glasses. The energy is built from a set of $M=\alpha N$ random functions
$V_{k}:\mathbf{S}^{N-1}\to\mathbb{R}$ that are centered Gaussians with
covariance

$\overline{V_{i}(\mathbf{x})V_{j}(\mathbf{x}^{\prime})}=\delta_{ij}f\left(\frac{\mathbf{x}\cdot\mathbf{x}^{\prime}}{N}\right)$
(70)

Each of the $V_{k}$ is an independent spherical spin glass. The total energy
is minus the sum of squares of the $V_{k}$, or

$H(\mathbf{x})=-\frac{1}{2}\sum_{k=1}^{M}V_{k}(\mathbf{x})^{2}$ (71)

The landscape complexity and large deviations of the ground state for the
least-squares version of this problem were recently studied in a linear
context, with $f(q)=\sigma^{2}+aq$ [41, 42, 43, 44]. Some results on the
ground state of the general nonlinear problem can also be found in [45], and a
solution to the equilibrium problem can be found in [46]. Those works indicate
that the low-lying minima of the least squares problem tend to be either
replica symmetric or full replica symmetry breaking. To avoid either a trivial
analysis or a very complex one, we instead focus on maximizing the sum of
squares, or minimizing (71).

The minima of (71) have a more amenable structure for study than the maxima,
as they are typically described by a 1rsb-like structure. There is a heuristic
intuition for this: in the limit of $M\to 1$, this problem is just minus 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 maxima are transformed into
minima. Therefore, the bottom of the landscape doesn’t qualitatively change.
The top, however, consists of the zero-energy level set in the spherical spin
glass. This level set is well-connected, and so the highest states should also
be well connected and flat.

Focusing on the bottom of the landscape and therefore dealing with a 1rsb-like
problem makes our analysis easier. Algorithms will tend to be stuck in the
ways they are in hard optimization problems, and we will be able to predict
where. Therefore, we will study the most squares problem rather than the least
squares one. We calculate the complexity of minima of (71) in Appendix D,
which corresponds to maximizing the sum of squares, under a replica symmetric
ansatz (which covers 1rsb-like problems) for arbitrary covariance $f$, and we
calculate the complexity of marginal minima in this section.

As in the previous sections, we used the method of Lagrange multipliers to
analyse stationary points on the constrained configuration space. The
Lagrangian and its associated gradient and Hessian are

$L(\mathbf{x},\omega)=-\frac{1}{2}\bigg{(}\sum_{k}^{M}V_{k}(\mathbf{x})^{2}-\omega\big{(}\|\mathbf{x}\|^{2}-N\big{)}\bigg{)}\\\
$ (72) $\displaystyle\nabla
H(\mathbf{x},\omega)=-\sum_{k}^{M}V_{k}(\mathbf{x})\partial
V_{k}(\mathbf{x})+\omega\mathbf{x}$ (73) $\displaystyle\begin{aligned}
&\operatorname{Hess}H(\mathbf{x},\omega)\\\
&\qquad=-\sum_{k}^{M}\left[\partial V_{k}(\mathbf{x})\partial
V_{k}(\mathbf{x})-V_{k}(\mathbf{x})\partial\partial
V_{k}(\mathbf{x})\right]+\omega I\end{aligned}$ (74)

As in the spherical and multispherical spin glasses, fixing the trace of the
Hessian is equivalent to constraining the value of the Lagrange multiplier
$\omega=\mu$.

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. Following the framework outlined in
Section III.3, the replicated number of stationary points conditioned on
energy $E$, trace $\mu$, and minimum eigenvalue $\lambda^{*}$ is given by

$\displaystyle\mathcal{N}(E,\mu,\lambda^{*})^{n}$ $\displaystyle=\int
d\hat{\beta}\,d\hat{\lambda}\prod_{a=1}^{n}\lim_{m_{a}\to
0}\prod_{\alpha=1}^{m_{a}}d\boldsymbol{\phi}_{a}^{\alpha}\exp\left\\{\delta^{\alpha
1}N(\hat{\beta}E+\hat{\lambda}\lambda^{*})-\frac{1}{2}\int
d1\,d2\,\left[B^{\alpha}(1,2)\sum_{k=1}^{M}V_{k}(\boldsymbol{\phi}_{a}^{\alpha}(1,2))^{2}-\mu\|\boldsymbol{\phi}_{a}^{\alpha}(1,2)\|^{2}\right]\right\\}$
(75)

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

$\delta\big{(}V_{k}(\boldsymbol{\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}(\boldsymbol{\phi}_{a}^{\alpha}(1,2))-v_{ka}^{\alpha}(1,2)\big{)}\right]$
(76)

where we have introduced auxiliary superfields $\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{v}V$ 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

$\overline{\exp\left[i\sum_{k}^{M}\sum_{a}^{n}\sum_{\alpha}^{m_{a}}\int
d1\,d2\,\hat{v}_{ka}^{\alpha}(1,2)V_{k}(\boldsymbol{\phi}_{a}^{\alpha}(1,2))\right]}=-\frac{1}{2}\sum_{ab}^{n}\sum_{\alpha\gamma}^{m_{a}}\sum_{k}^{M}\int
d1\,d2\,d3\,d4\,\hat{v}_{ka}^{\alpha}(1,2)f\big{(}\boldsymbol{\phi}_{a}^{\alpha}(1,2)\cdot\boldsymbol{\phi}_{b}^{\gamma}(3,4)\big{)}\hat{v}_{kb}^{\gamma}(3,4)$
(77)

The entire integrand is now factorized in the indices $k$ and quadratic in the
superfields $v$ and $\hat{v}$ with the kernel

$\begin{bmatrix}B^{\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{(}\boldsymbol{\phi}_{a}^{\alpha}(1,2)\cdot\boldsymbol{\phi}_{b}^{\gamma}(3,4)\big{)}\end{bmatrix}$
(78)

The integration over $v$ and $\hat{v}$ results in a term in the effective
action of the form

$\displaystyle-\frac{M}{2}\log\operatorname{sdet}\bigg{[}\delta(1,3)\,\delta(2,4)\delta_{ab}\delta^{\alpha\gamma}$
(79) $\displaystyle\hskip
70.0001pt+B^{\alpha}(1,2)f\big{(}\boldsymbol{\phi}_{a}^{\alpha}(1,2)\cdot\boldsymbol{\phi}_{b}^{\gamma}(3,4)\big{)}\bigg{]}$

When expanded, the supermatrix
$\boldsymbol{\phi}_{a}^{\alpha}(1,2)\cdot\boldsymbol{\phi}_{b}^{\gamma}(3,4)$
is constructed of the scalar products of the real and Grassmann vectors that
make up $\boldsymbol{\phi}$. The change of variables to these order parameters
again results in the Jacobian of (49), contributing

$\frac{N}{2}\log\det J-\frac{N}{2}\log\det G^{2}$ (80)

to the effective action.

Up to this point, the expressions are general and independent of a given
ansatz. However, we expect that the order parameters $X$ and $\hat{X}$ are
zero, since again we are in a setting with no signal or external field.
Applying this ansatz here avoids a dramatically more complicated expression
for the effective action. 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 this ansatz, taking the superdeterminant in (79) yields

$\displaystyle-\frac{M}{2}\log\det\left\\{\left[f^{\prime}(C)\odot
D-\hat{\beta}I+\left(R^{\circ 2}-G^{\circ
2}+I\sum_{\alpha\gamma}2(\delta^{\alpha 1}\hat{\lambda}+\beta)(\delta^{\gamma
1}\hat{\lambda}+\beta)(Q^{\alpha\gamma})^{2}\right)\odot
f^{\prime\prime}(C)\right]f(C)+(I-R\odot f^{\prime}(C))^{2}\right\\}$
$\displaystyle\hskip
160.00024pt-n\frac{M}{2}\log\det\big{[}\delta_{\alpha\gamma}-2(\delta_{\alpha
1}\hat{\lambda}+\beta)Q^{\alpha\gamma}\big{]}+M\log\det\big{[}I+G\odot
f^{\prime}(C)\big{]}$ (81)

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 into its own term. This is the realization of the fact
that the Hessian is 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
functions. This gives the ansatz

$\displaystyle C=I$ $\displaystyle R=rI$ $\displaystyle D=dI$ $\displaystyle
G=gI$ (82)

We further take a planted replica symmetric structure for the matrix $Q$,
identical to that in (9). This results in

$\displaystyle\Sigma_{\lambda^{*}}(E,\mu)=\frac{1}{N}\lim_{n\to
0}\frac{\partial}{\partial n}\int
d\hat{\beta}\,d\hat{\lambda}\,dr\,dd\,dg\,dq_{0}\,d\tilde{q}_{0}\,$ (83)
$\displaystyle\hskip 80.00012pt\times
e^{nN\mathcal{S}_{\mathrm{RSS}}(\hat{\beta},\hat{\lambda},r,d,g,q_{0},\tilde{q}_{0}\mid\lambda^{*},E,\mu,\beta)}$

with an effective action

$\displaystyle\mathcal{S}_{\mathrm{RSS}}(\hat{\beta},\hat{\lambda},r,d,g,q_{0},\tilde{q}_{0}\mid\lambda^{*},E,\mu,\beta)=\hat{\beta}E-\mu(r+g+\hat{\lambda})+\hat{\lambda}\lambda^{*}+\frac{1}{2}\log\left(\frac{d+r^{2}}{g^{2}}\times\frac{1-2q_{0}+\tilde{q}_{0}^{2}}{(1-q_{0})^{2}}\right)$
(84)
$\displaystyle\quad-\frac{\alpha}{2}\log\Bigg{(}\frac{1-4f^{\prime}(1)\big{[}\beta(1-q_{0})+\frac{1}{2}\hat{\lambda}-\beta(\beta+\hat{\lambda})(1-2q_{0}+\tilde{q}_{0}^{2})f^{\prime}(1)\big{]}}{\big{[}1-2(1-q_{0})\beta
f^{\prime}(1)\big{]}^{2}}$
$\displaystyle\qquad\qquad\qquad\times\frac{f(1)\big{[}f^{\prime}(1)d-\hat{\beta}-f^{\prime\prime}(1)\big{(}r^{2}-g^{2}+4q_{0}^{2}\beta^{2}-4\tilde{q}_{0}^{2}\beta(\beta+\hat{\lambda})+4\beta\hat{\lambda}+2\hat{\lambda}^{2}\big{)}\big{]}+(1-rf^{\prime}(1))^{2}}{\big{[}1+gf^{\prime}(1)\big{]}^{2}}\Bigg{)}$

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 like in (15)
and (16). Upon making this substitution and taking the zero-temperature limit,
we find

$\displaystyle\mathcal{S}_{\mathrm{RSS}}(\hat{\beta},\hat{\lambda},r,d,g,y,\Delta
z\mid\lambda^{*},E,\mu,\infty)=\hat{\beta}E-\mu(r+g+\hat{\lambda})+\hat{\lambda}\lambda^{*}+\frac{1}{2}\log\left(\frac{d+r^{2}}{g^{2}}\times\frac{y^{2}-2\Delta
z}{y^{2}}\right)$ (85)
$\displaystyle-\frac{\alpha}{2}\log\left(\frac{1-2(2y+\hat{\lambda})f^{\prime}(1)+4(y^{2}-2\Delta
z)f^{\prime}(1)^{2}}{\big{[}1-2yf^{\prime}(1)\big{]}^{2}}\times\frac{f(1)\big{[}f^{\prime}(1)d-\hat{\beta}-f^{\prime\prime}(1)(r^{2}-g^{2}+8(y\hat{\lambda}+\Delta
z)+2\hat{\lambda}^{2})\big{]}+\big{[}1-rf^{\prime}(1)\big{]}^{2}}{\big{[}1+gf^{\prime}(1)\big{]}^{2}}\right)$

We can finally write the complexity with fixed energy $E$, stability $\mu$,
and minimum eigenvalue $\lambda^{*}$ as

$\displaystyle\Sigma_{\lambda^{*}}(E,\mu)$ (86)
$\displaystyle=\underset{\hat{\beta},\hat{\lambda},r,d,g,y,\Delta
z}{\operatorname{extremum}}\mathcal{S}_{\textrm{RSS}}(\hat{\beta},\hat{\lambda},r,d,g,y,\Delta
z\mid\lambda^{*},E,\mu,\infty)$

Note that unlike the previous two examples, the effective action in this case
does not split into two largely independent pieces, one relating to the
eigenvalue problem and one relating to the ordinary complexity. Instead, the
order parameters related to the eigenvalue problem are mixed throughout the
effective action with those of the ordinary complexity. This is a signal of
the fact that the sum of squares problem is not Gaussian, while the previous
two examples are. In all non-Gaussian problems, conditioning on properties of
the Hessian cannot be done independently from the complexity, and the method
introduced in this paper becomes necessary.

Figure 3:  Dominant and marginal complexity in the nonlinear sum of squares
problem for $\alpha=\frac{3}{2}$ 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.

The marginal complexity can be derived from (86) using the condition (32) to
fix $\mu$ to the marginal stability $\mu_{\textrm{m}}(E)$ and then evaluating
the complexity at that stability as in (33). Fig. 3 shows the marginal
complexity in a sum-of-squares model with $\alpha=\frac{3}{2}$ and
$f(q)=q^{2}+q^{3}$. Also shown is the dominant complexity computed in Appendix
D. As the figure demonstrates, the range of energies at which marginal minima
are found can differ significantly from those implied by the dominant
complexity, with the lowest energy significantly higher than the ground state
and the highest energy significantly higher than the threshold.

Figure 4:  The stability, or shift of the trace, for dominant and marginal
optima in the nonlinear sum of squares problem for $\alpha=\frac{3}{2}$ and
$f(q)=q^{2}+q^{3}$.

Fig. 4 shows the associated marginal stability $\mu_{\mathrm{m}}(E)$ for the
same model. Recall that the definition of the marginal stability in (32) is
that which eliminates the variation of $\Sigma_{\lambda^{*}}(E,\mu)$ with
respect to $\lambda^{*}$ at the point $\lambda^{*}=0$. Unlike in the Gaussian
spherical spin glass, in this model $\mu_{\mathrm{m}}(E)$ varies with energy
in a nontrivial way. The figure also shows the dominant stability, which is
the stability associated with the dominant complexity and coincides with the
marginal stability only at the threshold energy.

Figure 5:  Real and imaginary parts of the complexity $\Sigma_{0}(E,\mu)$ with
fixed minimum eigenvalue $\lambda^{*}=0$ as a function of $\mu$ in the
nonlinear sum of squares problem with $\alpha=\frac{3}{2}$,
$f(q)=q^{2}+q^{3}$, and $E\simeq-6.47$. The vertical line depicts the value of
the marginal stability $\mu_{\mathrm{m}}$.

Because this version of the model has no signal, we were able to use the
heuristic (32) to fix the marginal stability. However, we could also have used
the more general method for finding a pseudogapped Hessian spectrum by
locating the value of $\mu$ at which the complexity develops an imaginary
part, as described in Section II.3 and pictured in Fig. 1. The real and
imaginary parts of the complexity $\Sigma_{0}(E,\mu)$ are plotted in Fig. 5 as
a function of $\mu$ at fixed energy. The figure also shows the marginal
stability $\mu_{\mathrm{m}}$ predicted by the variational approach (32). The
marginal stability corresponds to precisely the point at which an imaginary
part develops in the complexity. This demonstrates that the principles we used
to determine the marginal stability continue to hold even in non-Gaussian
cases where the complexity and the condition to fix the minimum eigenvalue are
tangled together.

In our companion paper, we use a sum of squared random functions model to
explore the relationship between the marginal complexity and the performance
of two generic algorithms: gradient descent and approximate message passing
[21]. We show that the range of energies where the marginal complexity is
positive does effectively bound the performance of these algorithms. At the
moment the comparison is restricted to models with small polynomial powers
appearing in $f(q)$ and with small $\alpha$ for computational reasons.
However, using the dmft results already found for these models it should be
possible to make comparisons in a wider family of models [47, 48].

The results for the marginal complexity are complimentary to rigorous results
on the performance of algorithms in the least squares case, which focus on
bounds for $\alpha$ and the parameters of $f$ necessary for zero-energy
solutions to exist and be found by algorithms [49, 50]. After more work to
evaluate the marginal complexity in the full rsb case, it will be interesting
to compare the bounds implied by the distribution of marginal minima with
those made by other means.

##  V Conclusions

We have introduced a method for conditioning complexity on the marginality of
stationary points. This method is general, and permits conditioning without
first needing to understand the statistics of the Hessian at stationary
points. We used our approach to study 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 marginal complexity in the third model of sums of
squared random functions can be used to effectively bound algorithmic
performance [21].

There are some limitations to the approach we relied on in this paper. The
main limitation is our restriction to signalless landscapes, where there is no
symmetry-breaking favored direction. This allowed us to treat stationary
points with isolated eigenvalues as atypical, and therefore find the marginal
stability $\mu_{\mathrm{m}}$ 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, 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
approach. It is 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.

Finally, the problem of predicting which marginal minima are able to attract
some dynamics and which cannot attract any dynamics looms large over this
work. As we discussed briefly at the end of Section IV.2, in some simple
contexts it is easy to see why certain marginal minima are not viable, but at
the moment we do not know how to generalize this. Ideas related to the self-
similarity and stochastic stability of minima have recently been suggested as
a route to understanding this problem, but this approach is still in its
infancy [51].

###### Acknowledgements.

JK-D is supported by a DynSysMath Specific Initiative of the INFN.

##  Appendix A A primer on superspace

In this appendix we review the algebra of superspace [52]. 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

$\int d\theta\,\theta=1\qquad\int d\theta\,1=0$ (87)

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 superspace function, then the integral of
$f$ with respect to a Grassmann index can be evaluated using this property of
the series expansion by

$\int d\theta\,f(a+b\theta)=\int
d\theta\,\left[f(a)+f^{\prime}(a)b\theta\right]=f^{\prime}(a)b$ (88)

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

$\boldsymbol{\phi}(1)=\mathbf{x}+\bar{\theta}_{1}\boldsymbol{\eta}+\bar{\boldsymbol{\eta}}\theta_{1}+\bar{\theta}_{1}\theta_{1}i\hat{\mathbf{x}}$
(89)

where $\mathbf{x},\hat{\mathbf{x}}\in\mathbb{R}^{N}$ and
$\bar{\boldsymbol{\eta}},\boldsymbol{\eta}$ are $N$-dimensional Grassmann
vectors. The dependence of $\boldsymbol{\phi}$ on 1 indicates the index of
Grassmann variables $\bar{\theta}_{1},\theta_{1}$ inside, since we will
sometimes want to use, e.g., $\boldsymbol{\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:

$\displaystyle f\big{(}\boldsymbol{\phi}(1)\big{)}$
$\displaystyle=f(\mathbf{x})+\big{(}\bar{\theta}_{1}\boldsymbol{\eta}+\bar{\boldsymbol{\eta}}\theta_{1}+\bar{\theta}_{1}\theta_{1}i\hat{\mathbf{x}}\big{)}^{T}\partial
f(\mathbf{x})$ (90)
$\displaystyle\quad+\frac{1}{2}\big{(}\bar{\theta}_{1}\boldsymbol{\eta}+\bar{\boldsymbol{\eta}}\theta_{1}\big{)}^{T}\partial\partial
f(\mathbf{x})\big{(}\bar{\theta}_{1}\boldsymbol{\eta}+\bar{\boldsymbol{\eta}}\theta_{1}\big{)}$
$\displaystyle=f(\mathbf{x})+\big{(}\bar{\theta}_{1}\boldsymbol{\eta}+\bar{\boldsymbol{\eta}}\theta_{1}+\bar{\theta}_{1}\theta_{1}i\hat{\mathbf{x}}\big{)}^{T}\partial
f(\mathbf{x})$
$\displaystyle\qquad-\bar{\theta}_{1}\theta_{1}\bar{\boldsymbol{\eta}}^{T}\partial\partial
f(\mathbf{x})\boldsymbol{\eta}$

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

$\int
d\theta_{1}\,d\bar{\theta}_{1}\,f\big{(}\boldsymbol{\phi}(1)\big{)}=i\hat{\mathbf{x}}^{T}\partial
f(\mathbf{x})-\bar{\boldsymbol{\eta}}^{T}\partial\partial
f(\mathbf{x})\boldsymbol{\eta}$ (91)

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

$\displaystyle\int d\mathbf{x}\,\delta\big{(}\nabla
H(\mathbf{x})\big{)}\,\det\operatorname{Hess}H(\mathbf{x})$ (92)
$\displaystyle\qquad=\int
d\mathbf{x}\,d\bar{\boldsymbol{\eta}}\,d\boldsymbol{\eta}\,\frac{d\hat{\mathbf{x}}}{(2\pi)^{N}}\,e^{i\hat{\mathbf{x}}^{T}\nabla
H(\mathbf{x})-\bar{\boldsymbol{\eta}}^{T}\operatorname{Hess}H(\mathbf{x})\boldsymbol{\eta}}$
$\displaystyle\qquad=\int d\boldsymbol{\phi}\,e^{\int
d1\,H(\boldsymbol{\phi}(1))}$

where we have written the measures $d1=d\theta_{1}\,d\bar{\theta}_{1}$ and
$d\boldsymbol{\phi}=d\mathbf{x}\,d\bar{\boldsymbol{\eta}}\,d\boldsymbol{\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

$\int d\boldsymbol{\phi}\,d\hat{\beta}\,e^{\hat{\beta}E+\int
d1\,(1-\hat{\beta}\bar{\theta}_{1}\theta_{1})H(\boldsymbol{\phi}(1))}$ (93)

which a small calculation confirms results in the same expression as (35).

The reason why this transformation is a simplification 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
$\boldsymbol{\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.,

$(M\boldsymbol{\phi})(1)=\int d1\,M(1,2)\boldsymbol{\phi}(2)$ (94)

The identity supermatrix is given by

$\delta(1,2)=(\bar{\theta}_{1}-\bar{\theta}_{2})(\theta_{1}-\theta_{2})I$ (95)

Integrals involving superfields contracted into such operators result in
schematically familiar expressions, like that of the standard Gaussian:

$\int
d\boldsymbol{\phi}\,e^{-\frac{1}{2}\int\,d1\,d2\,\boldsymbol{\phi}(1)^{T}M(1,2)\boldsymbol{\phi}(2)}=(\operatorname{sdet}M)^{-1/2}$
(96)

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

$\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}$ (97)

where each of the blocks is a $2N\times 2N$ real matrix. Then the
superdeterminant of $M$ is given by

$\operatorname{sdet}M=\det(A-BD^{-1}C)\det(D)^{-1}$ (98)

which is the same as the normal expression for the determinant of a block
matrix save for the inverse of $\det D$. Likewise, the supertrace of $M$ is is
given by

$\operatorname{sTr}M=\operatorname{Tr}A-\operatorname{Tr}D$ (99)

The same method can be used to calculate the superdeterminant and supertrace
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}\\}$.

##  Appendix B BRST symmetry

When the trace $\mu$ is not fixed, there is an unusual symmetry in the
dominant complexity of minima [30, 29, 31]. This arises from considering the
Kac–Rice formula as a kind of gauge fixing procedure [53]. Around each
stationary point consider making the coordinate transformation
$\mathbf{u}=\nabla H(\mathbf{x})$. Then, in the absence of fixing the trace of
the Hessian to $\mu$, the Kac–Rice measure becomes

$\int d\nu(\mathbf{x},\boldsymbol{\omega}\mid
E)=\int\sum_{\sigma}d\mathbf{u}\,\delta(\mathbf{u})\,\delta\big{(}NE-H(\mathbf{x}_{\sigma})\big{)}$
(100)

where the sum is over stationary points $\sigma$. 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 function has been dropped and Grassmann vectors
introduced to represent the determinant of the Hessian, this symmetry can be
simplified considerably. Due to the expansion properties of Grassmann
integrals, any appearance of $-\bar{\boldsymbol{\eta}}\boldsymbol{\eta}^{T}$
in the integrand resolves to $(\operatorname{Hess}H)^{-1}$. The symmetry of
the measure can then be written

$\mathbf{x}\mapsto\mathbf{x}-\bar{\boldsymbol{\eta}}\boldsymbol{\eta}^{T}\delta\mathbf{u}=\mathbf{x}+\bar{\boldsymbol{\eta}}\delta\epsilon$
(101)

where $\delta\epsilon=-\boldsymbol{\eta}^{T}\delta\mathbf{u}$ is a Grassmann
number. This establishes that
$\delta\mathbf{x}=\bar{\boldsymbol{\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$. This gives

$\displaystyle\delta\mathbf{x}=\bar{\boldsymbol{\eta}}\,\delta\epsilon$
$\displaystyle\delta\hat{\mathbf{x}}=-i\hat{\beta}\bar{\boldsymbol{\eta}}\,\delta\epsilon$
$\displaystyle\delta\boldsymbol{\eta}=-i\hat{\mathbf{x}}\,\delta\epsilon$
$\displaystyle\delta\bar{\boldsymbol{\eta}}=0$ (102)

so that the differential form of the symmetry is

$\mathcal{D}=\bar{\boldsymbol{\eta}}\cdot\frac{\partial}{\partial\mathbf{x}}-i\hat{\beta}\bar{\boldsymbol{\eta}}\cdot\frac{\partial}{\partial\hat{\mathbf{x}}}-i\hat{\mathbf{x}}\cdot\frac{\partial}{\partial\boldsymbol{\eta}}$
(103)

The Ward identities associated with this symmetry give rise to relationships
among the order parameters. These identities come from applying the
differential symmetry to Grassmann-valued order parameters, and are

$\displaystyle\begin{aligned}
0&=\frac{1}{N}\mathcal{D}\langle\mathbf{x}_{a}\cdot\boldsymbol{\eta}_{b}\rangle=\frac{1}{N}\left[\langle\bar{\boldsymbol{\eta}}_{a}\cdot\boldsymbol{\eta}_{b}\rangle-i\langle\mathbf{x}_{a}\cdot\hat{\mathbf{x}}_{b}\rangle\right]\\\
&=G_{ab}+R_{ab}\end{aligned}$ (104) $\displaystyle\begin{aligned}
0&=\frac{i}{N}\mathcal{D}\langle\hat{\mathbf{x}}_{a}\cdot\boldsymbol{\eta}_{b}\rangle=\frac{1}{N}\left[\hat{\beta}\langle\bar{\boldsymbol{\eta}}_{a}\cdot\boldsymbol{\eta}_{b}\rangle+\langle\hat{\mathbf{x}}_{a}\cdot\hat{\mathbf{x}}_{b}\rangle\right]\\\
&=\hat{\beta}G_{ab}+D_{ab}\end{aligned}$ (105)

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.

##  Appendix C Spectral density in the multispherical spin glass

In this appendix we derive an expression for the asymptotic spectral density
of the Hessian in the two-sphere multispherical spin glass that we describe in
Section IV.2. We use a typical approach of employing replicas to compute the
resolvent [54]. The resolvent for the Hessian of the multispherical model is
given by an integral over
$\mathbf{y}=[\mathbf{y}^{(1)},\mathbf{y}^{(2)}]\in\mathbb{R}^{2N}$ as

$\displaystyle G(\lambda)$ $\displaystyle=\lim_{n\to
0}\int\|\mathbf{y}_{1}\|^{2}\,\prod_{a=1}^{n}d\mathbf{y}_{a}\,\exp\left\\{-\frac{1}{2}\mathbf{y}_{a}^{T}(\operatorname{Hess}H(\mathbf{x},\boldsymbol{\omega})-\lambda
I)\mathbf{y}_{a}\right\\}$ (106) $\displaystyle=\lim_{n\to
0}\int\big{(}\|\mathbf{y}_{1}^{(1)}\|^{2}+\|\mathbf{y}_{1}^{(2)}\|^{2}\big{)}\,\prod_{a=1}^{n}d\mathbf{y}_{a}\,\exp\left\\{-\frac{1}{2}\begin{bmatrix}\mathbf{y}_{a}^{(1)}\\\
\mathbf{y}_{a}^{(2)}\end{bmatrix}^{T}\left(\begin{bmatrix}\partial_{1}\partial_{1}H_{1}(\mathbf{x}^{(1)})+\omega_{1}I&-\epsilon
I\\\ -\epsilon
I&\partial_{2}\partial_{2}H_{2}(\mathbf{x}^{(2)})+\omega_{2}I\end{bmatrix}-\lambda
I\right)\begin{bmatrix}\mathbf{y}_{a}^{(1)}\\\
\mathbf{y}_{a}^{(2)}\end{bmatrix}\right\\}$

If $Y_{ab}^{(ij)}=\frac{1}{N}\mathbf{y}_{a}^{(i)}\cdot\mathbf{y}_{b}^{(j)}$ is
the matrix of overlaps of the vectors $\mathbf{y}$, then a short and standard
calculation involving the average over $H$ and the change of variables from
$\mathbf{y}$ to $Y$ yields

$\overline{G(\lambda)}=N\lim_{n\to 0}\int
dY\,\big{(}Y_{11}^{(11)}+Y_{11}^{(22)}\big{)}\,e^{nN\mathcal{S}(Y)}$ (107)

where the effective action $\mathcal{S}$ is given by

$\displaystyle\mathcal{S}(Y)=\lim_{n\to
0}\frac{1}{n}\Bigg{\\{}\frac{1}{4}\sum_{ab}^{n}\left[f_{1}^{\prime\prime}(1)(Y_{ab}^{(11)})^{2}+f_{2}^{\prime\prime}(1)(Y_{ab}^{(22)})^{2}\right]+\frac{1}{2}\sum_{a}^{n}\left[2\epsilon
Y_{aa}^{(12)}+(\lambda-\omega_{1})Y_{aa}^{(11)}+(\lambda-\omega_{2})Y_{aa}^{(22)}\right]$
(108) $\displaystyle+\frac{1}{2}\log\det\begin{bmatrix}Y^{(11)}&Y^{(12)}\\\
Y^{(12)}&Y^{(22)}\end{bmatrix}\Bigg{\\}}$

Making the replica symmetric ansatz $Y_{ab}^{(ij)}=y^{(ij)}\delta_{ab}$ for
each of the matrices $Y^{(ij)}$ yields

$\displaystyle\mathcal{S}(y)$
$\displaystyle=\frac{1}{4}\left[f_{1}^{\prime\prime}(1)(y^{(11)})^{2}+f_{2}^{\prime\prime}(1)(y^{(22)})^{2}\right]+\epsilon
y^{(12)}$ (109)
$\displaystyle\qquad+\frac{1}{2}\left[(\lambda-\omega_{1})y^{(11)}+(\lambda-\omega_{2})y^{(22)}\right]$
$\displaystyle\qquad+\frac{1}{2}\log(y^{(11)}y^{(22)}-y^{(12)}y^{(12)})$

while the average resolvent becomes

$\overline{G(\lambda)}=N(y^{(11)}+y^{(22)})$ (110)

for $y^{(11)}$ and $y^{(22)}$ evaluated at a saddle point of $\mathcal{S}$.
The spectral density at large $N$ is then given by the discontinuity in its
imaginary point on the real axis, or

$\rho(\lambda)=\frac{1}{2\pi
iN}\left(\overline{G(\lambda+i0^{+})}-\overline{G(\lambda+i0^{-})}\right)$
(111)

##  Appendix D Complexity of dominant optima for sums of squared random
functions

Here we share an outline of the derivation of formulas for the complexity of
dominant optima in sums of squared random functions of section IV.3. While in
this paper we only treat problems with a replica symmetric structure, formulas
for the effective action are generic to any rsb structure and provide a
starting point for analyzing the challenging full rsb setting.

Using the $\mathbb{R}^{N|2}$ superfields

$\boldsymbol{\phi}_{a}(1)=\mathbf{x}_{a}+\bar{\theta}_{1}\boldsymbol{\eta}_{a}+\bar{\boldsymbol{\eta}}_{a}\theta_{1}+\bar{\theta}_{1}\theta_{1}\hat{\mathbf{x}}_{a},$
(112)

the replicated count of stationary points can be written

$\displaystyle\mathcal{N}(E,\mu)^{n}=\int
d\hat{\beta}\prod_{a=1}^{n}\,d\boldsymbol{\phi}_{a}\,\exp\bigg{[}N\hat{\beta}E$
(113) $\displaystyle\qquad-\frac{1}{2}\int
d1\,\left(B(1)\sum_{k=1}^{M}V_{k}(\boldsymbol{\phi}_{a}(1))^{2}-\mu\|\boldsymbol{\phi}_{a}(1)\|^{2}\right)\bigg{]}$

for $B(1)=1-\hat{\beta}\bar{\theta}_{1}\theta_{1}$. The derivation of the
complexity follows from here nearly identically to that in Appendix A.2 of
Fyodorov and Tublin with superoperations replacing standard ones [43]. First
we insert Dirac $\delta$ functions to fix each of the $M$ energies
$V_{k}(\boldsymbol{\phi}_{a}(1))$ as

$\displaystyle\delta\big{(}V_{k}(\boldsymbol{\phi}_{a}(1))-v_{ka}(1)\big{)}$
(114) $\displaystyle=\int d\hat{v}_{ka}\,\exp\left[i\int
d1\,\hat{v}_{ka}(1)\big{(}V_{k}(\boldsymbol{\phi}_{a}(1))-v_{ka}(1)\big{)}\right]$

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 (114), which is linear. The average over the disorder
can then be computed, which yields

$\displaystyle\overline{\exp\left[i\sum_{k=1}^{M}\sum_{a=1}^{n}\int
d1\,\hat{v}_{ka}(1)V_{k}(\boldsymbol{\phi}_{a}(1))\right]}$ (115)
$\displaystyle=\exp\left[-\frac{1}{2}\sum_{k=1}^{M}\sum_{ab}^{n}\int
d1\,d2\,\hat{v}_{ka}(1)f\left(\frac{\boldsymbol{\phi}_{a}(1)\cdot\boldsymbol{\phi}_{b}(2)}{N}\right)\hat{v}_{kb}(2)\right]$

The result is factorized in the indices $k$ and Gaussian in the superfields
$v$ and $\hat{v}$ with kernel

$\begin{bmatrix}B(1)\delta_{ab}\delta(1,2)&i\delta_{ab}\delta(1,2)\\\
i\delta_{ab}\delta(1,2)&f\left(\frac{\boldsymbol{\phi}_{a}(1)\cdot\boldsymbol{\phi}_{b}(2)}{N}\right)\end{bmatrix}$
(116)

Making the $M$ independent Gaussian integrals, we find

$\displaystyle\mathcal{N}(E,\mu)^{n}=\int
d\hat{\beta}\left(\prod_{a=1}^{n}d\boldsymbol{\phi}_{a}\right)\exp\bigg{[}nN\hat{\beta}E+\frac{\mu}{2}\sum_{a}^{n}\int
d1\,\|\boldsymbol{\phi}_{a}\|^{2}$ (117)
$\displaystyle\quad-\frac{M}{2}\log\operatorname{sdet}\left(\delta_{ab}\delta(1,2)+B(1)f\left(\frac{\boldsymbol{\phi}_{a}(1)\cdot\boldsymbol{\phi}_{b}(2)}{N}\right)\right)\bigg{]}$

We make a change of variables from the fields $\boldsymbol{\phi}$ to matrices
$\mathbb{Q}_{ab}(1,2)=\frac{1}{N}\boldsymbol{\phi}_{a}(1)\cdot\boldsymbol{\phi}_{b}(2)$.
This transformation results in a change of measure of the form

$\prod_{a=1}^{n}d\boldsymbol{\phi}_{a}=d\mathbb{Q}\,(\operatorname{sdet}\mathbb{Q})^{\frac{N}{2}}=d\mathbb{Q}\,\exp\left[\frac{N}{2}\log\operatorname{sdet}\mathbb{Q}\right]$
(118)

We therefore have

$\displaystyle\mathcal{N}(E,\mu)^{n}=\int
d\hat{\beta}\,d\mathbb{Q}\,\exp\bigg{\\{}nN\hat{\beta}E+N\frac{\mu}{2}\operatorname{sTr}\mathbb{Q}+\frac{N}{2}\log\operatorname{sdet}\mathbb{Q}-\frac{M}{2}\log\operatorname{sdet}\left[\delta_{ab}\delta(1,2)+B(1)f(\mathbb{Q}_{ab}(1,2))\right]\bigg{\\}}$
(119)

We now need to blow up our supermatrices into our physical order parameters.
We have from the definition of $\boldsymbol{\phi}$ and $\mathbb{Q}$ that

$\displaystyle\mathbb{Q}_{ab}(1,2)=C_{ab}-G_{ab}(\bar{\theta}_{1}\theta_{2}+\bar{\theta}_{2}\theta_{1})-R_{ab}(\bar{\theta}_{1}\theta_{1}+\bar{\theta}_{2}\theta_{2})-D_{ab}\bar{\theta}_{1}\theta_{2}\bar{\theta}_{2}\theta_{2}$
(120)

where $C$, $R$, $D$, and $G$ are the matrices defined in (47). Other possible
combinations involving scalar products between fermionic and bosonic variables
do not contribute at physical saddle points [33]. Inserting this expansion
into the expression above and evaluating the superdeterminants and supertrace,
we find

$\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)}$
(121)

where the effective action is given by

$\displaystyle\mathcal{S}_{\mathrm{KR}}(\hat{\beta},C,R,D,G)$
$\displaystyle=\hat{\beta}E+\lim_{n\to
0}\frac{1}{n}\Bigg{(}-\mu\operatorname{Tr}(G+R)+\frac{1}{2}\log\det\big{[}G^{-2}(CD+R^{2})\big{]}+\alpha\log\det\big{[}I+G\odot
f^{\prime}(C)\big{]}$ (122)
$\displaystyle\qquad-\frac{\alpha}{2}\log\det\left[\Big{(}f^{\prime}(C)\odot
D-\hat{\beta}I+(G^{\circ 2}-R^{\circ 2})\odot
f^{\prime\prime}(C)\Big{)}f(C)+(I-R\odot f^{\prime}(C))^{2}\right]\Bigg{)}$

where $\odot$ gives the Hadamard or componentwise product between the matrices
and $A^{\circ n}$ gives the Hadamard power of $A$, while other products and
powers are matrix products and powers.

In the case where $\mu$ is not specified, we can make use of the BRST symmetry
of Appendix B whose Ward identities give $D=\hat{\beta}R$ and $G=-R$. Using
these relations, the effective action becomes particularly simple:

$\mathcal{S}_{\mathrm{KR}}(\hat{\beta},C,R)=\hat{\beta}E+\frac{1}{2}\lim_{n\to
0}\frac{1}{n}\Big{(}\log\det(I+\hat{\beta}CR^{-1})-\alpha\log\det\left[I-\hat{\beta}f(C)\big{(}I-R\odot
f^{\prime}(C)\big{)}^{-1}\right]\Big{)}$ (123)

This effective action is general for arbitrary matrices $C$ and $R$, and
therefore arbitrary rsb order. 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

$\displaystyle\mathcal{S}_{\mathrm{KR}}(\hat{\beta},c_{0},r,r_{0})=\hat{\beta}E+\frac{1}{2}\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]$
(124)
$\displaystyle\qquad-\frac{\alpha}{2}\left[\log\left(1-\frac{\hat{\beta}\big{(}f(1)-f(c_{0})\big{)}}{1-rf^{\prime}(1)+r_{0}f^{\prime}(c_{0})}\right)-\frac{\hat{\beta}f(c_{0})+r_{0}f^{\prime}(c_{0})}{1-\hat{\beta}\big{(}f(1)-f(c_{0})\big{)}-rf^{\prime}(1)+rf^{\prime}(c_{0})}+\frac{r_{0}f^{\prime}(c_{0})}{1-rf^{\prime}(1)+r_{0}f^{\prime}(c_{0})}\right]$

When $f(0)=0$ as in the cases directly studied in this work, this further
simplifies as $c_{0}=r_{0}=0$. The effective action is then

$\mathcal{S}_{\mathrm{KR}}(\hat{\beta},r)=\hat{\beta}E+\frac{1}{2}\log\left(1+\frac{\hat{\beta}}{r}\right)-\frac{\alpha}{2}\log\left(1-\frac{\hat{\beta}f(1)}{1-rf^{\prime}(1)}\right)$
(125)

Extremizing this expression with respect to the order parameters $\hat{\beta}$
and $r$ produces the red line of dominant minima shown in Fig. 3.

## References

  * Ros and Fyodorov [2023] V. Ros and Y. V. Fyodorov, The high-dimensional landscape paradigm: Spin-glasses, and beyond, in  _Spin Glass Theory and Far Beyond_ (World Scientific, 2023) pp. 95–114.
  * Biroli and Garrahan [2013] G. Biroli and J. P. Garrahan, Perspective: The glass transition, The Journal of Chemical Physics 138, 12A301 (2013).
  * Krzakala and Kurchan [2007] F. Krzakala and J. Kurchan, Landscape analysis of constraint satisfaction problems, Physical Review E 76, 021122 (2007).
  * Folena  _et al._ [2020] G. Folena, S. Franz, and F. Ricci-Tersenghi, Rethinking mean-field glassy dynamics and its relation with the energy landscape: The surprising case of the spherical mixed $p$-spin model, Physical Review X 10, 031045 (2020).
  * Folena and Zamponi [2023] G. Folena and F. Zamponi, On weak ergodicity breaking in mean-field spin glasses, SciPost Physics 15, 109 (2023).
  * El Alaoui and Montanari [2020] A. El Alaoui and A. Montanari, Algorithmic thresholds in mean field spin glasses, arXiv:2009.11481v1 [cond-mat.stat-mech] (2020).
  * Kent-Dobias and Kurchan [2023] J. Kent-Dobias and J. Kurchan, How to count in hierarchical landscapes: a full solution to mean-field complexity, Physical Review E 107, 064111 (2023).
  * Müller and Wyart [2015] M. Müller and M. Wyart, Marginal stability in structural, spin, and electron glasses, Annual Review of Condensed Matter Physics 6, 177 (2015).
  * Anderson [1984] P. W. Anderson, Lectures on amorphous systems, in  _Ill-Condensed Matter_ (North-Holland Publishing Company, 1984) pp. 159–261.
  * Sommers and Dupont [1984] H. J. Sommers and W. Dupont, Distribution of frozen fields in the mean-field theory of spin glasses, Journal of Physics C: Solid State Physics 17, 5785 (1984).
  * Parisi [1995] G. Parisi, On the statistical properties of the large time zero temperature dynamics of the SK model, arXiv:cond-mat/9501045v1 (1995).
  * Horner [2007] H. Horner, Time dependent local field distribution and metastable states in the SK-spin-glass, The European Physical Journal B 60, 413 (2007).
  * Pankov [2006] S. Pankov, Low-temperature solution of the Sherrington-Kirkpatrick model, Physical Review Letters 96, 197204 (2006).
  * Erba  _et al._ [2024] V. Erba, F. Behrens, F. Krzakala, and L. Zdeborová, Quenches in the Sherrington-Kirkpatrick model, arXiv:2405.04267v2 (2024).
  * Efros and Shklovskii [1985] A. L. Efros and B. I. Shklovskii, Coulomb interaction in disordered systems with localized electronic states, in  _Electron–Electron Interactions in Disordered Systems_, Modern Problems in Condensed Matter Sciences, Vol. 10 (Elsevier, 1985) pp. 409–482.
  * Shklovskii [2024] B. I. Shklovskii, Half century of Efros-Shklovskii Coulomb gap. romance with Coulomb interaction and disorder, arXiv:2403.19793v5 (2024).
  * Folena and Urbani [2022] G. Folena and P. Urbani, Marginal stability of soft anharmonic mean field spin glasses, Journal of Statistical Mechanics: Theory and Experiment 2022, 053301 (2022).
  * Monasson [1995] R. Monasson, Structural glass transition and the entropy of the metastable states, Physical Review Letters 75, 2847 (1995).
  * Fyodorov [2004] Y. V. Fyodorov, Complexity of random energy landscapes, glass transition, and absolute value of the spectral determinant of random matrices, Physical Review Letters 92, 240601 (2004).
  * Bray and Dean [2007] A. J. Bray and D. S. Dean, Statistics of critical points of Gaussian fields on large-dimensional spaces, Physical Review Letters 98, 150201 (2007).
  * Kent-Dobias [2024a] J. Kent-Dobias, Algorithm-independent bounds on complex optimization through the statistics of marginal optima, arXiv:2407.02092 [cond-mat.dis-nn] (2024a).
  * Ikeda [2023] H. Ikeda, Bose–Einstein-like condensation of deformed random matrix: a replica approach, Journal of Statistical Mechanics: Theory and Experiment 2023, 023302 (2023).
  * Kent-Dobias [2024b] J. Kent-Dobias, Arrangement of nearby minima and saddles in the mixed spherical energy landscapes, SciPost Physics 16, 001 (2024b).
  * Fyodorov and Le Doussal [2013] Y. V. Fyodorov and P. Le Doussal, Topology trivialization and large deviations for the minimum in the simplest random optimization, Journal of Statistical Physics 154, 466 (2013).
  * Fyodorov and Le Doussal [2018] Y. V. Fyodorov and P. Le Doussal, Hessian spectrum at the global minimum of high-dimensional random landscapes, Journal of Physics A: Mathematical and Theoretical 51, 474002 (2018).
  * Dean and Majumdar [2006] D. S. Dean and S. N. Majumdar, Large deviations of extreme eigenvalues of random matrices, Physical Review Letters 97, 160201 (2006).
  * Kac [1943] M. Kac, On the average number of real roots of a random algebraic equation, Bulletin of the American Mathematical Society 49, 314 (1943).
  * Rice [1939] S. O. Rice, The distribution of the maxima of a random curve, American Journal of Mathematics 61, 409 (1939).
  * Annibale  _et al._ [2003a] A. Annibale, A. Cavagna, I. Giardina, and G. Parisi, Supersymmetric complexity in the Sherrington-Kirkpatrick model, Physical Review E 68, 061103 (2003a).
  * Annibale  _et al._ [2003b] A. Annibale, A. Cavagna, I. Giardina, G. Parisi, and E. Trevigne, The role of the Becchi–Rouet–Stora–Tyutin supersymmetry in the calculation of the complexity for the Sherrington–Kirkpatrick model, Journal of Physics A: Mathematical and General 36, 10937 (2003b).
  * Annibale  _et al._ [2004] A. Annibale, G. Gualdi, and A. Cavagna, Coexistence of supersymmetric and supersymmetry-breaking states in spherical spin-glasses, Journal of Physics A: Mathematical and General 37, 11311 (2004).
  * Crisanti  _et al._ [1993] A. Crisanti, H. Horner, and H.-J. Sommers, The spherical $p$-spin interaction spin-glass model, Zeitschrift für Physik B Condensed Matter 92, 257 (1993).
  * Kurchan [1992] J. Kurchan, Supersymmetry in spin glass dynamics, Journal de Physique I 2, 1333 (1992).
  * Subag [2021] E. Subag, TAP approach for multi-species spherical spin glasses I: general theory, arXiv:2111.07132v1 (2021).
  * Subag [2023] E. Subag, TAP approach for multispecies spherical spin glasses II: The free energy of the pure models, The Annals of Probability 51, 1004 (2023).
  * Bates and Sohn [2022a] E. Bates and Y. Sohn, Crisanti–Sommers formula and simultaneous symmetry breaking in multi-species spherical spin glasses, Communications in Mathematical Physics 394, 1101 (2022a).
  * Bates and Sohn [2022b] E. Bates and Y. Sohn, Free energy in multi-species mixed $p$-spin spherical models, Electronic Journal of Probability 27, 1 (2022b).
  * Huang and Sellke [2023a] B. Huang and M. Sellke, Strong topological trivialization of multi-species spherical spin glasses, arXiv:2308.09677v2 (2023a).
  * Huang and Sellke [2023b] B. Huang and M. Sellke, Algorithmic threshold for multi-species spherical spin glasses, arXiv:2303.12172v2 (2023b).
  * Huang and Sellke [2024] B. Huang and M. Sellke, Optimization algorithms for multi-species spherical spin glasses, Journal of Statistical Physics 191, 29 (2024).
  * Fyodorov [2019] Y. V. Fyodorov, A spin glass model for reconstructing nonlinearly encrypted signals corrupted by noise, Journal of Statistical Physics 175, 789 (2019).
  * Fyodorov and Tublin [2020] Y. V. Fyodorov and R. Tublin, Counting stationary points of the loss function in the simplest constrained least-square optimization, Acta Physica Polonica B 51, 1663 (2020).
  * Fyodorov and Tublin [2022] Y. V. Fyodorov and R. Tublin, Optimization landscape in the simplest constrained random least-square problem, Journal of Physics A: Mathematical and Theoretical 55, 244008 (2022).
  * Vivo [2024] P. Vivo, Random linear systems with quadratic constraints: from random matrix theory to replicas and back, arXiv:2401.03209v2 (2024).
  * Tublin [2022] R. Tublin,  _A Few Results in Random Matrix Theory and Random Optimization_ , Ph.D. thesis, King’s College London (2022).
  * Urbani [2023] P. Urbani, A continuous constraint satisfaction problem for the rigidity transition in confluent tissues, Journal of Physics A: Mathematical and Theoretical 56, 115003 (2023).
  * Kamali and Urbani [2023a] P. J. Kamali and P. Urbani, Dynamical mean field theory for models of confluent tissues and beyond, SciPost Physics 15, 219 (2023a).
  * Kamali and Urbani [2023b] P. J. Kamali and P. Urbani, Stochastic gradient descent outperforms gradient descent in recovering a high-dimensional signal in a glassy energy landscape, arXiv:2309.04788v2 (2023b).
  * Montanari and Subag [2023] A. Montanari and E. Subag, Solving overparametrized systems of random equations: I. model and algorithms for approximate solutions, arXiv:2306.13326v1 (2023).
  * Montanari and Subag [2024] A. Montanari and E. Subag, On Smale’s 17th problem over the reals, arXiv:2405.01735v1 (2024).
  * Urbani [2024] P. Urbani, Statistical physics of complex systems: glasses, spin glasses, continuous constraint satisfaction problems, high-dimensional inference and neural networks, arXiv:2405.06384v1 (2024).
  * DeWitt [1992] B. S. DeWitt,  _Supermanifolds_ , 2nd ed., Cambridge monographs on mathematical physics (Cambridge University Press, Cambridge ; New York, 1992).
  * Zinn-Justin [2002] J. Zinn-Justin,  _Quantum field theory and critical phenomena_ (Clarendon Press Oxford University Press, Oxford New York, 2002).
  * Livan  _et al._ [2018] G. Livan, M. Novaes, and P. Vivo,  _Introduction to Random Matrices_, SpringerBriefs in Mathematical Physics, Vol. 26 (Springer International Publishing, 2018).