\documentclass[submission, Phys]{SciPost}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,latexsym,graphicx}
\usepackage[bitstream-charter]{mathdesign}
\usepackage[dvipsnames]{xcolor}
\usepackage{anyfontsize,authblk}

\urlstyle{sf}

% Fix \cal and \mathcal characters look (so it's not the same as \mathscr)
\DeclareSymbolFont{usualmathcal}{OMS}{cmsy}{m}{n}
\DeclareSymbolFontAlphabet{\mathcal}{usualmathcal}

\hypersetup{
  colorlinks=true,
  urlcolor={blue!50!black},
  citecolor={blue!50!black},
  filecolor={blue!50!black},
  linkcolor={blue!50!black}
}

\title{
  On the topology of solutions to random continuous constraint satisfaction problems
}

\author{Jaron Kent-Dobias\footnote{\url{jaron.kent-dobias@roma1.infn.it}}}
\affil{Istituto Nazionale di Fisica Nucleare, Sezione di Roma I, Italy}

\begin{document}

\maketitle

\begin{abstract}
  We consider the set of solutions to $M$ random polynomial equations with
  independent Gaussian coefficients and a target value $V_0$ on the $(N-1)$-sphere. When solutions
  exist, they form a manifold. We compute the average Euler characteristic of
  this manifold in the limit of large $N$, and find different behavior
  depending on the ratio $\alpha=M/N$. When $\alpha<\alpha_\text{onset}$, the
  average characteristic is 2 and there is a single connected component, while
  for $\alpha_\text{onset}<\alpha<\alpha_\text{shatter}$ a large connected
  component coexists with many disconnected components. When
  $\alpha>\alpha_\text{shatter}$ the large connected component vanishes, and the
  entire manifold vanishes for $\alpha>\alpha_\text{\textsc{sat}}$. In the
  limit $\alpha\to0$ there is a correspondence between this problem and the
  topology of constant-energy level sets in the spherical spin glasses. We
  conjecture that the energy $E_\text{shatter}$ associated with the vanishing of
  the large connected component corresponds to the asymptotic limit of gradient
  descent from a random initial condition.
\end{abstract}

\tableofcontents

\section{Introduction}

Constraint satisfaction problems seek configurations that simultaneously
satisfy a set of equations, and form a basis for thinking about problems as
diverse as neural networks \cite{Mezard_2009_Constraint}, granular materials
\cite{Franz_2017_Universality}, ecosystems \cite{Altieri_2019_Constraint}, and
confluent tissues \cite{Urbani_2023_A}. All but the last of these examples deal
with sets of inequalities, while the last considers a set of equality
constraints. Inequality constraints are familiar in situations like zero-cost
solutions in neural networks with ReLu activations and stable equilibrium in the
forces between physical objects. Equality constraints naturally appear in the
zero-gradient solutions to overparameterized smooth neural networks and in vertex models of tissues.

In such problems, there is great interest in characterizing structure in the
set of solutions, which can be influential in how algorithms behave when trying
to solve them \cite{Baldassi_2016_Unreasonable, Baldassi_2019_Properties,
Beneventano_2023_On}. Here, we show how \emph{topological} information about
the set of solutions can be calculated in a simple model of satisfying random
nonlinear equalities. This allows us to reason about the connectivity of this
solution set. The topological properties revealed by this calculation yield
surprising results for the well-studied spherical spin glasses, where a
topological transition thought to occur at a threshold energy $E_\text{th}$
where marginal minima are dominant is shown to occur at a different energy
$E_\text{shatter}$. We conjecture that this difference resolves an outstanding
problem in gradient descent dynamics in these systems.

We consider the problem of finding configurations $\mathbf x\in\mathbb R^N$
lying on the $(N-1)$-sphere $\|\mathbf x\|^2=N$ that simultaneously satisfy $M$
nonlinear constraints $V_k(\mathbf x)=V_0$ for $1\leq k\leq M$ and some
constant $V_0\in\mathbb R$. The nonlinear constraints are taken to be centered
Gaussian random functions with covariance
\begin{equation} \label{eq:covariance}
  \overline{V_i(\mathbf x)V_j(\mathbf x')}
  =\delta_{ij}f\left(\frac{\mathbf x\cdot\mathbf x'}N\right)
\end{equation}
for some choice of function $f$. When the covariance function $f$ is polynomial, the
$V_k$ are also polynomial, with a term of degree $p$ in $f$ corresponding to
all possible terms of degree $p$ in $V_k$. In particular, taking
\begin{equation}
  V_k(\mathbf x)
  =\sum_{p=0}^\infty\frac1{p!}\sqrt{\frac{f^{(p)}(0)}{N^p}}
  \sum_{i_1\cdots i_p}^NJ^{(k,p)}_{i_1\cdots i_p}x_{i_1}\cdots x_{i_p}
\end{equation}
with the elements of the tensors $J^{(k,p)}$ as independently distributed
unit normal random variables satisfies \eqref{eq:covariance}. The size of the
series coefficients of $f$ therefore control the variances in the coefficients
of random polynomial constraints. When $M=1$, this problem corresponds to the
level set of a spherical spin glass with energy density $E=\sqrt{N}V_0$.


This problem or small variations thereof have attracted attention recently for
their resemblance to encryption, optimization, and vertex models of confluent
tissues \cite{Fyodorov_2019_A, Fyodorov_2020_Counting,
  Fyodorov_2022_Optimization, Urbani_2023_A, Kamali_2023_Dynamical,
  Kamali_2023_Stochastic, Urbani_2024_Statistical, Montanari_2023_Solving,
Montanari_2024_On, Kent-Dobias_2024_Conditioning, Kent-Dobias_2024_Algorithm-independent}. In each of these cases, the authors studied properties of
the cost function
\begin{equation} \label{eq:cost}
  \mathscr C(\mathbf x)=\frac12\sum_{k=1}^M\big[V_k(\mathbf x)-V_0\big]^2
\end{equation}
which achieves zero only for configurations that satisfy all the constraints.
Here we dispense with the cost function and study the set of solutions
directly.
This set
can be written as
\begin{equation}
  \Omega=\big\{\mathbf x\in\mathbb R^N\mid \|\mathbf x\|^2=N,V_k(\mathbf x)=V_0
  \;\forall\;k=1,\ldots,M\big\}
\end{equation}
Because the constraints are all smooth functions, $\Omega$ is almost always a manifold without singular points.\footnote{The conditions for a singular point are that
$0=\frac\partial{\partial\mathbf x}V_k(\mathbf x)$ for all $k$. This is
equivalent to asking that the constraints $V_k$ all have a stationary point at
the same place. When the $V_k$ are independent and random, this is vanishingly
unlikely, requiring $NM+1$ independent equations to be simultaneously satisfied.
This means that different connected components of the set of solutions do not
intersect, nor are there self-intersections, without extraordinary fine-tuning.}
We study the topology of the manifold $\Omega$ by two related means: its
average Euler characteristic, and the average number of stationary points of a
linear height function restricted to the manifold. These measures tell us
complementary pieces of information, respectively the alternating sum and
direct sum of the Betti numbers of $\Omega$. We find that for the varied cases
we study, these two always coincide at the largest exponential order in $N$,
putting strong constraints on the resulting topology and geometry.

\section{Methods}

\subsection{The average Euler characteristic}

The Euler characteristic $\chi$ of a manifold is a topological invariant \cite{Hatcher_2002_Algebraic}. It is
perhaps most familiar in the context of connected compact orientable surfaces, where it
characterizes the number of handles in the surface: $\chi=2(1-\#)$ for $\#$
handles. For general $d$, the Euler characteristic of the $d$-sphere is $2$ if $d$ is even and 0 if $d$ is odd. The canonical method for computing the Euler characteristic is done by
defining a complex on the manifold in question, essentially a
higher-dimensional generalization of a polygonal tiling. Then $\chi$ is given
by an alternating sum over the number of cells of increasing dimension, which
for 2-manifolds corresponds to the number of vertices, minus the number of
edges, plus the number of faces.

Morse theory offers another way to compute the Euler characteristic using the
statistics of stationary points of a function $H:\Omega\to\mathbb R$ \cite{Audin_2014_Morse}. For
functions $H$ without any symmetries with respect to the manifold, the surfaces
of gradient flow between adjacent stationary points form a complex. The
alternating sum over cells to compute $\chi$ becomes an alternating sum over
the count of stationary points of $H$ with increasing index, or
\begin{equation}
  \chi=\sum_{i=0}^N(-1)^i\mathcal N_H(\text{index}=i)
\end{equation}
Conveniently, we can express this abstract sum as an integral over the manifold
using a small variation on the Kac--Rice formula for counting stationary
points. Since the sign of the determinant of the Hessian matrix of $H$ at a
stationary point is equal to its index, if we count stationary points including
the sign of the determinant, we arrive at the Euler characteristic, or
\begin{equation} \label{eq:kac-rice}
  \chi(\Omega)=\int_\Omega d\mathbf x\,\delta\big(\nabla H(\mathbf x)\big)\det\operatorname{Hess}H(\mathbf x)
\end{equation}
When the Kac--Rice formula is used to \emph{count} stationary points, the sign
of the determinant is a nuisance that one must take pains to preserve
\cite{Fyodorov_2004_Complexity}. Here we are correct to exclude it.

We need to choose a function $H$ for our calculation. Because $\chi$ is
a topological invariant, any choice will work so long as it does not share some
symmetry with the underlying manifold, i.e., that it $H$ satisfies the Smale condition. Because our manifold of random
constraints has no symmetries, we can take a simple height function $H(\mathbf
x)=\mathbf x_0\cdot\mathbf x$ for some $\mathbf x_0\in\mathbb R^N$ with
$\|\mathbf x_0\|^2=N$. $H$ is a height function because when $\mathbf x_0$ is
used as the polar axis, $H$ gives the height on the sphere.

We treat the integral over the implicitly defined manifold $\Omega$ using the
method of Lagrange multipliers. We introduce one multiplier $\omega_0$ to
enforce the spherical constraint and $M$ multipliers $\omega_k$ to enforce the vanishing of
each of the $V_k$, resulting in the Lagrangian
\begin{equation} \label{eq:lagrangian}
  L(\mathbf x,\pmb\omega)
  =H(\mathbf x)+\frac12\omega_0\big(\|\mathbf x\|^2-N\big)
  +\sum_{k=1}^M\omega_k\big(V_k(\mathbf x)-V_0\big)
\end{equation}
The integral over $\Omega$ in \eqref{eq:kac-rice} then becomes
\begin{equation} \label{eq:kac-rice.lagrange}
  \chi(\Omega)=\int_{\mathbb R^N} d\mathbf x\int_{\mathbb R^{M+1}}d\pmb\omega
  \,\delta\big(\partial L(\mathbf x,\pmb\omega)\big)
  \det\partial\partial L(\mathbf x,\pmb\omega)
\end{equation}
where $\partial=[\frac\partial{\partial\mathbf x},\frac\partial{\partial\pmb\omega}]$
is the vector of partial derivatives with respect to all $N+M+1$ variables.
This integral is now in a form where standard techniques from mean-field theory
can be applied to calculate it.

In order for certain Gaussian integrals in the following calculation to be
well-defined, it is necessary to treat instead the Lagrangian problem above
with $\pmb\omega\mapsto i\pmb\omega$. This transformation does not effect the
Dirac $\delta$ functions of the gradient, but it does change the determinant by
a factor of $i^{N+M+1}$. We will see that the result of the rest of the
calculation neglecting this factor is real. Since the Euler characteristic is
also necessarily real, this indicates an inconsistency with this transformation
when $N+M+1$ is odd. In fact, the Euler characteristic is always zero for
odd-dimensional manifolds. This is the signature of it in this problem.

\subsubsection{Calculation of the average Euler characteristic}

To evaluate the average of $\chi$ over the constraints, we first translate the $\delta$ functions and determinant to integral form, with
\begin{align}
  \delta\big(\partial L(\mathbf x,\pmb\omega)\big)
  &=\int\frac{d\hat{\mathbf x}}{(2\pi)^N}\frac{d\hat{\pmb\omega}}{(2\pi)^{M+1}}
  e^{i[\hat{\mathbf x},\hat{\pmb\omega}]\cdot\partial L(\mathbf x,\pmb\omega)}
  \\
  \det\partial\partial L(\mathbf x,\pmb\omega)
  &=\int d\bar{\pmb\eta}\,d\pmb\eta\,d\bar{\pmb\gamma}\,d\pmb\gamma\,
  e^{-[\bar{\pmb\eta},\bar{\pmb\gamma}]^T\partial\partial L(\mathbf x,\pmb\omega)[\pmb\eta,\pmb\gamma]}
\end{align}
where $\hat{\mathbf x}$ and $\hat{\pmb\omega}$ are ordinary vectors and
$\bar{\pmb\eta}$, $\pmb\eta$, $\bar{\pmb\gamma}$, and $\pmb\gamma$ are
Grassmann vectors. With these expressions substituted into
\eqref{eq:kac-rice.lagrange}, the result is a integral over an exponential
whose argument is linear in the random functions $V_k$. These functions can
therefore be averaged over, and the resulting expression treated with standard
methods. Details of this calculation can be found in Appendix~\ref{sec:euler}.
The result is the reduction of the average Euler characteristic to an expression of the
form
\begin{equation}
  \overline{\chi(\Omega)}
  =\int dR\,dD\,dm\,d\hat m\,g(R,D,m,\hat m)\,e^{N\mathcal S_\Omega(R,D,m,\hat m)}
\end{equation}
where $g$ is a prefactor subexponential in $N$, and $\mathcal S_\Omega$ is an effective action defined by
\begin{equation}
  \begin{aligned}
    \mathcal S_\Omega(R,D,m,\hat m\mid\alpha,V_0)
    &=\hat m-\frac\alpha2\left[
      \log\left(1+\frac{f(1)D}{f'(1)R^2}\right)
      +\frac{V_0^2}{f(1)}\left(1+\frac{f'(1)R^2}{f(1)D}\right)^{-1}
    \right] \\
    &\hspace{7em}+\frac12\log\left(
      1+\frac{(1-m^2)D+\hat m^2-2Rm\hat m}{R^2}
    \right)
  \end{aligned}
\end{equation}
The remaining order parameters defined by the scalar products
\begin{align}
  R=-i\frac1N\mathbf x\cdot\hat{\mathbf x}
  &&
  D=\frac1N\hat{\mathbf x}\cdot\hat{\mathbf x}
  &&
  m=\frac1N\mathbf x\cdot\mathbf x_0
  &&
  \hat m=-i\frac1N\hat {\mathbf x}\cdot\mathbf x_0
\end{align}

This integral can be evaluated by a saddle point method. For reasons we will
see, it is best to extremize with respect to $R$, $D$, and $\hat m$, leaving a
new effective action of $m$ alone. This can be solved to give
\begin{equation}
  D=-\frac{m+R^*(m)}{1-m^2} \qquad \hat m=0
\end{equation}
\begin{equation}
  \begin{aligned}
    R^*(m)
    &=\frac{-m(1-m^2)}{2[f(1)-(1-m^2)f'(1)]^2}
    \Bigg[
      \alpha V_0^2f'(1)
      +(2-\alpha)f(1)\left(\frac{f(1)}{1-m^2}-f'(1)\right) \\
    &\quad+\operatorname{sgn}(m)\alpha\sqrt{
      \tfrac{4V_0^2}\alpha f(1)f'(1)\left[\tfrac{f(1)}{1-m^2}-f'(1)\right]
      +\left[\tfrac{f(1)^2}{1-m^2}-\big(V_0^2+f(1)\big)f'(1)\right]^2
    }
    \Bigg]
  \end{aligned}
\end{equation}
\begin{equation}
  \mathcal S_\Omega(m)
  =-\frac\alpha2\bigg[
    \log\left(
      1-\frac{f(1)}{f'(1)}\frac{1+m/R^*}{1-m^2}
    \right)
    +\frac{V_0^2}{f(1)}\left(
      1-\frac{f'(1)}{f(1)}\frac{1-m^2}{1+m/R^*}
    \right)^{-1}
  \bigg]
  +\frac12\log\left(-\frac m{R^*}\right)
\end{equation}
To finish evaluating the integral, this expression should be maximized with
respect to $m$. The order parameter $m$ is both physical and interpretable, as
it gives the overlap of the configuration $\mathbf x$ with the height axis
$\mathbf x_0$. Therefore, the value $m^*$ that maximizes this action can be
understood as the latitude on the sphere where most of the contribution to the
Euler characteristic is made.

The action $\mathcal S_\Omega$ is extremized with respect to $m$ at $m^*=0$ or $m^*=-R^*$.
In the latter case, $m^*$ takes the value
\begin{equation}
  m^*=\pm\sqrt{1-\frac{\alpha}{f'(1)}\big(V_0^2+f(1)\big)}
\end{equation}
and $\mathcal S_\Omega(m^*)=0$. However, when
\begin{equation}
  V_0^2>V_\text{on}^2\equiv\frac{1-\alpha+\sqrt{1-\alpha}}\alpha f(1)
\end{equation}
$R^*(m^*)$ becomes complex and the solution is no longer valid. Likewise, when
\begin{equation}
  V_0^2<V_\text{sh}^2\equiv\frac{2(1+\sqrt{1-\alpha})-\alpha}{\alpha}f(1)\left(1-\frac{f(1)}{f'(1)}\right)
\end{equation}
the maximum at $m=0$ becomes complex and that solution is invalid. Examples of
$\mathcal S_\Omega$ in difference regimes are shown in Fig.~\ref{fig:action}.
These transition values of the target $V_0$ correspond with transition values in $\alpha$ of
\begin{align}
  \alpha_\text{on}=1-\left(\frac{V_0^2}{V_0^2+f(1)}\right)^2
  &&
  \alpha_\text{sh}=4V_0^2f(1)f'(1)\frac{f'(1)-f(1)}{\big((V_0^2+f(1))f'(1)-f(1)^2\big)^2}
\end{align}

\begin{figure}
  \includegraphics{figs/action_1.pdf}
  \hspace{-3.5em}
  \includegraphics{figs/action_3.pdf}

  \caption{
    The effective action governing the average Euler characteristic as a function of the overlap
    $m=\frac1N\mathbf x\cdot\mathbf x_0$ with the height direction for two
    different homogeneous polynomial functions and a variety of target values $V_0$. In both
    plots $\alpha=\frac12$. \textbf{Left:} With linear functions there are two
    regimes. For small $V_0$, there are maxima at $m=\pm m^*$ where the action
    is zero, while after the satisfiability transition at $V_0=V_\text{\textsc{sat}}=1$, $m^*$
    goes to zero and the action becomes negative. \textbf{Right:} With nonlinear
    functions there are four regimes. For small $V_0$ the behavior is the same
    as in the linear case, with zero action. After an onset transition at
    $V_0=V_\text{on}\simeq1.099$ the maxima are at the edge of validity of the
    action and the action is positive. At a shattering transition at
    $V_0=V_\text{sh}\simeq1.394$, $m^*$ goes to zero and the action is positive.
    Finally, at the satisfiability transition
    $V_0=V_\text{\textsc{sat}}\simeq1.440$ the action becomes negative.
  } \label{fig:action}
\end{figure}




\subsubsection{How to interpret the average Euler characteristic}

It is not straightforward to directly use the average Euler characteristic to
infer something about the number of connected components in the set of
solutions. To understand why, a simple example is helpful. Consider the set of
solutions on the sphere $\|\mathbf x\|^2=N$ that satisfy the single quadratic
constraint
\begin{equation}
  0=\sum_{i=1}^N\sigma_ix_i^2
\end{equation}
where each $\sigma_i$ is taken to be $\pm1$ with equal probability. If we take $\mathbf x$ to be ordered such that all terms with $\sigma_i=+1$ come first, this gives
\begin{equation}
  0=\sum_{i=1}^{N_+}x_i^2-\sum_{i=N_++1}^Nx_i^2
\end{equation}
where $N_+$ is the number of terms with $\sigma_i=+1$. The topology of the resulting manifold can be found by adding and subtracting this constraint from the spherical one, which gives
\begin{align}
  \frac12=\sum_{i=1}^{N_+}x_i^2
  \qquad
  \frac12=\sum_{i=N_++1}^{N}x_i^2
\end{align}
These are two independent equations for spheres of radius $1/\sqrt2$, one of
dimension $N_+$ and the other of dimension $N-N_+$. Therefore, the topology of
the configuration space is that of $S^{N_+-1}\times S^{N-N_+-1}$. The Euler
characteristic of a product space is the product of the Euler characteristics,
and so we have $\chi(\Omega)=\chi(S^{N_+-1})\chi(S^{N-N_+-1})$.

What is the average value of the Euler characteristic over values of
$\sigma_i$? First, recall that the Euler characteristic of a sphere $S^d$ is 2
when $d$ is even and 0 when $d$ is odd. When $N$ is odd, any value
of $N_+$ will result in one of the two spheres in the product to be
odd-dimensional, and therefore $\chi(\Omega)=0$, as is always true for
odd-dimensional manifolds. When $N$ is even, there are two possibilities: when $N_+$ is even then both spheres are odd-dimensional, while when $N_+$ is odd then both spheres are even-dimensional.
The number of terms $N_+$ with $\sigma_i=+1$ is distributed with the binomial distribution
\begin{equation}
  P(N_+)=\frac1{2^N}\binom{N}{N_+}
\end{equation}
Therefore, the average Euler characteristic for even $N$ is
\begin{equation}
  \overline{\chi(\Omega)}
  =\sum_{N_+=0}^NP(N_+)\chi(S^{N_+-1})\chi(S^{N-N_+-1})
  =\frac4{2^N}\sum_{n=0}^{N/2}\binom{N}{2n}
  =2
\end{equation}
Thus we find the average Euler characteristic in this simple example is 2
despite the fact that the possible manifolds resulting from the constraints
have characteristics of either 0 or 4.

The many positive values the Euler characteristic can take is one problem, but
a larger one is the many negative values it can take. A manifold with many
handles will have a very negative Euler characteristic, and our calculation
relying on the saddle point of an exponential integral would be invalid. In the
regime of $\mathcal S_\Omega(m)$ where the action is complex-valued, is this
breakdown the result of a negative average characteristic, or is it the result
of another reason? It is difficult to answer this with the previous calculation
alone. However, in the next subsection we will make a complementary calculation
that rules out the negative Euler characteristic picture. Instead, we will see
that in the complex region, the large-deviation function in $m$ that $\mathcal
S_\Omega(m)$ represents breaks down due to the vanishingly small probability of
finding any stationary points.

\subsection{Complexity of a linear test function}

One way to definitely narrow possible interpretations of the average Euler
characteristic is to compute a complementary average. The Euler characteristic
is the alternating sum of numbers of critical points of different index. If
instead we make the direct sum
\begin{equation}
  \mathcal N_H(\Omega)=\sum_{i=0}^N\mathcal N_H(\text{index}=i)
  =\int_\Omega d\mathbf x\,\delta\big(\nabla H(\mathbf x)\big)
  \,\big|\det\operatorname{Hess}H(\mathbf x)\big|
\end{equation}
we find the total number of stationary points. The formula is exactly the same
as that for the average Euler characteristic except for an absolute value sign
around the determinant of the Hessian.

Understanding the number of stationary points as a function of latitude $m$
will clarify the meaning of our effective action for the average Euler
characteristic. This is because the average number of stationary points is a
nonnegative number. If the region of complex $\mathcal S_\Omega$ has a
well-defined number of stationary points, it indicates that we are looking at a
situation with a negative average Euler characteristic. On the other hand, if
the average number of stationary points yields a complex value at some latitude
$m$, it must be because it is either to large or small in $N$ to be captured by
the calculation, e.g., that it behaves like $e^{N^2\Sigma}$ or
$e^{-N^2\Sigma}$. The following calculation indicates this second situation:
the region of complex action is due to a lack of stationary points to
contribute to the Euler characteristic at those latitudes.

To compute the complexity, we follow a similar procedure to the Euler
characteristic. The main difference lies in how we treat the absolute value
function around the determinant. Following \cite{Fyodorov_2004_Complexity}, we
make use of the identity
\begin{equation}
  \begin{aligned}
  |\det A|
  &=\lim_{\epsilon\to0}\frac{(\det A)^2}{\sqrt{\det(A+i\epsilon I)}\sqrt{\det(A-i\epsilon I)}}
  \\
  &=\frac1{(2\pi)^N}\int d\bar{\pmb\eta}_1\,d\pmb\eta_1\,d\bar{\pmb\eta}_2\,d\pmb\eta_2\,d\mathbf a\,d\mathbf b\,
  e^{
    -\bar{\pmb\eta}_1^TA\pmb\eta_1-\bar{\pmb\eta}_2^TA\pmb\eta_2
    -\frac12\mathbf a^T(A+i\epsilon I)\mathbf a-\frac12\mathbf b^T(A-i\epsilon I)\mathbf b
  }
  \end{aligned}
\end{equation}
for an $N\times N$ matrix $A$. Here $\bar{\pmb\eta}_1$, $\pmb\eta_1$,
$\bar{\pmb\eta}_2$, and $\pmb\eta_2$ are Grassmann vectors and $\mathbf a$ and
$\mathbf b$ are regular vectors. This introduces many new order parameters into
the problem, but this is a difficulty of scale rather than principle. With this
identity substituted for the usual determinant one, the problem can be solved
much as before. The details of this solution are relegated to
Appendix~\ref{sec:complexity.details}. The result is that, to largest order in
$N$, the logarithm of the average number of stationary points is the same as
the logarithm of the average Euler characteristic.

\section{Results}

\subsection{Topology of solutions to many equations and the satisfiability transition}

The results of the previous sections indicate the following picture for the
topology of the solution manifold.

\begin{figure}
  \includegraphics[width=0.245\textwidth]{figs/connected.pdf}
  \includegraphics[width=0.245\textwidth]{figs/coexist.pdf}
  \includegraphics[width=0.245\textwidth]{figs/shattered.pdf}
  \includegraphics[width=0.245\textwidth]{figs/gone.pdf}

  \includegraphics{figs/bar.pdf}

  \caption{
    Cartoon of the topology of the solution manifold implied by our
    calculation. The arrow shows the vector $\mathbf x_0$ defining the height
    function. For $V_0<V_\text{on}$, the manifold has a single connected
    component. Above the onset with $V_\text{on}<V_0<V_\text{sh}$, the manifold
    has a large connected component around the equator, and many disconnected
    pieces in a certain range of latitudes. Above the shattering transition, or
    $V_\text{sh}<V_0<V_\text{\textsc{sat}}$, the large connected
    component vanishes and small disconnected pieces occupy the entire equatorial
    region. Finally, above the satisfiability transition
    $V_\text{\textsc{sat}}$ the manifold vanishes.
  } \label{fig:cartoons}
\end{figure}

\begin{figure}
  \includegraphics{figs/phases_1.pdf}
  \hspace{-3em}
  \includegraphics{figs/phases_2.pdf}
  \hspace{-3em}
  \includegraphics{figs/phases_3.pdf}

  \caption{
    Topological phases of the model for three different homogeneous covariance
    functions. The onset transition $V_\text{onset}$, shattering transition
    $V_\text{shatter}$, and satisfiability transition $V_\text{\textsc{sat}}$
    are indicated when they exist. In the limit of $\alpha\to0$, the behavior
    of level sets of the spherical spin glasses are recovered: the final plot
    shows how in the pure cubic model the ground state energy $E_\text{gs}$ and threshold energy
    $E_\text{th}$ correspond with the limits of the satisfiability and
    shattering transitions, respectively. Note that for mixed models with
    inhomogeneous covariance functions, $E_\text{th}$ is not the lower limit of
    $V_\text{sh}$.
  }
\end{figure}

\subsection{Topology of level sets of the spherical spin glasses and the dynamic threshold}

As indicated earlier, for $M=1$ the solution manifold corresponds to the energy
level set of a spherical spin glass with energy density $E=\sqrt NV_0$. All the
results from the previous sections follow, and can be translated to the spin
glasses by taking the limit $\alpha\to0$ while scaling $V_0=\sqrt\alpha E$. With a little algebra this procedure yields
\begin{equation}
  E_\text{onset}=\pm\sqrt{2f(1)}
\end{equation}
\begin{equation}
  E_\text{shatter}=\pm\sqrt{4f(1)\left(1-\frac{f(1)}{f'(1)}\right)}
\end{equation}
for the energies at which level sets of the spherical spin glasses have
disconnected pieces appear, and that at which a large connected component
vanishes. For the pure $p$-spin spherical spin glasses with $f(q)=\frac12q^p$,
$E_\text{shatter}=\sqrt{2(p-1)/p}$, precisely the threshold energy in these
models. This is expected, since threshold energy, defined as the place where
marginal minima are dominant in the landscape, is widely understood as the
place where level sets are broken into pieces.

However, for general mixed models the threshold energy is
\begin{equation}
  E_\mathrm{th}=\pm\frac{f(1)[f''(1)-f'(1)]+f'(1)^2}{f'(1)\sqrt{f''(1)}}
\end{equation}
which generally satisfies $|E_\text{shatter}|\leq|E_\text{th}|$.

\begin{figure}
  \includegraphics{figs/dynamics_2.pdf}
  \hspace{-0.5em}
  \includegraphics{figs/dynamics_3.pdf}

  \caption{
    Comparison of the shattering energy $E_\text{sh}$ with the asymptotic
    performance of gradient descent from a random initial condition in $p+s$
    models with $p=2$ and $p=3$ and varying $s$. The values of $\lambda$ depend on $p$ and $s$ and are taken from \cite{Folena_2023_On}. The points show the asymptotic performance
    extrapolated using two different methods and have unknown uncertainty, from \cite{Folena_2023_On}. Also
    shown is the annealed threshold energy $E_\text{th}$, where marginal minima
    are the most common type of stationary point. The section of $E_\text{sh}$
    that is dashed on the left plot indicates the continuation of the annealed
    result, whereas the solid portion gives the value calculated with a
    {\oldstylenums 1}\textsc{frsb} ansatz.
  }
\end{figure}


\section{Conclusions}


\cite{Franz_2016_The, Franz_2017_Universality, Franz_2019_Critical, Annesi_2023_Star-shaped, Baldassi_2023_Typical}

\paragraph{Acknowledgements}
The authors thank Pierfrancesco Urbani for helpful conversations on these topics.

\paragraph{Funding information}
JK-D is supported by a \textsc{DynSysMath} Specific Initiative of the INFN.

\appendix

\section{Details of the calculation of the average Euler characteristic}
\label{sec:euler}

To make the calculation compact, we introduce superspace coordinates. Define the supervectors
\begin{equation}
  \pmb\phi(1)=\mathbf x+\bar\theta_1\pmb\eta+\bar{\pmb\eta}\theta_1+\bar\theta_1\theta_1i\hat{\mathbf x}
  \qquad
  \sigma_k(1)=\omega_k+\bar\theta_1\gamma_k+\bar\gamma_k\theta_1+\bar\theta_1\theta_1\hat\omega_k
\end{equation}
The Euler characteristic can be expressed using these supervectors as
\begin{equation}
  \begin{aligned}
    &\chi(\Omega)
    =\int d\pmb\phi\,d\pmb\sigma\,e^{\int d1\,L(\pmb\phi(1),\pmb\sigma(1))} \\
    &=\int d\pmb\phi\,d\pmb\sigma\,\exp\left\{
      \int d1\left[
        H\big(\pmb\phi(1)\big)
        +\frac i2\sigma_0(1)\left(\|\pmb\phi(1)\|^2-N\right)
        +i\sum_{k=1}^M\sigma_k(1)\Big(V_k\big(\pmb\phi(1)\big)-V_0\Big)
      \right]
    \right\}
  \end{aligned}
\end{equation}
Since this is an exponential integrand linear in the functions $V_k$, we can average over the functions to find
\begin{equation}
  \begin{aligned}
    \overline{\chi(\Omega)}
    =\int d\pmb\phi\,d\pmb\sigma\,\exp\Bigg\{
      \int d1\left[
        H(\pmb\phi(1))
        +\frac{i}2\sigma_0(1)\big(\|\pmb\phi(1)\|^2-N\big)
        -iV_0\sum_{k=1}^M\sigma_k(1)
      \right] \\
      -\frac12\int d1\,d2\,\sum_{k=1}^M\sigma_k(1)\sigma_k(2)f\left(\frac{\pmb\phi(1)\cdot\pmb\phi(2)}N\right)
    \Bigg\}
  \end{aligned}
\end{equation}
This is a Gaussian integral in the Lagrange multipliers with $1\leq k\leq M$.
Performing that integral yields
\begin{equation}
  \begin{aligned}
    \overline{\chi(\Omega)}
    &=\int d\pmb\phi\,d\sigma_0\,\exp\Bigg\{
      \int d1\left[
        H(\pmb\phi(1))
        +\frac{i}2\sigma_0(1)\big(\|\pmb\phi(1)\|^2-N\big)
      \right] \\
    &\hspace{5em}-\frac M2V_0^2\int d1\,d2\,f\left(\frac{\pmb\phi(1)\cdot\pmb\phi(2)}N\right)^{-1}
      -\frac M2\log\operatorname{sdet}f\left(\frac{\pmb\phi(1)\cdot\pmb\phi(2)}N\right)
    \Bigg\}
  \end{aligned}
\end{equation}
The supervector $\pmb\phi$ enters this expression as a function only of the
scalar product with itself and with the vector $\mathbf x_0$ inside the
function $H$. We therefore change variables to the superoperator $\mathbb Q$ and the supervector $\mathbb M$ defined by
\begin{equation}
  \mathbb Q(1,2)=\frac{\pmb\phi(1)\cdot\pmb\phi(2)}N
  \qquad
  \mathbb M(1)=\frac{\pmb\phi(1)\cdot\mathbf x_0}N
\end{equation}
These new variables can replace $\pmb\phi$ in the integral using a generalized Hubbard--Stratonovich transformation, which yields
\begin{equation}
  \begin{aligned}
    \overline{\chi(\Omega)}
    &=\int d\mathbb Q\,d\mathbb M\,d\sigma_0\,
    \left[g(\mathbb Q,\mathbb M)+O(N^{-1})\right]
    \,\exp\Bigg\{
      N\int d1\left[
        \mathbb M(1)
        +\frac{i}2\sigma_0(1)\big(\mathbb Q(1,1)-1\big)
      \right] \\
    &\hspace{5em}-\frac M2V_0^2\int d1\,d2\,f(\mathbb Q)^{-1}(1,2)
      -\frac M2\log\operatorname{sdet}f(\mathbb Q)
      +\frac N2\log\operatorname{sdet}(\mathbb Q-\mathbb M\mathbb M^T)
    \Bigg\}
  \end{aligned}
\end{equation}
where $g$ is a function of $\mathbb Q$ and $\mathbb M$ independent of $N$ and
$M$, detailed in Appendix~\ref{sec:prefactor}. To move on from this expression,
we need to expand the superspace notation. We can write
\begin{equation}
  \begin{aligned}
    \mathbb Q(1,2)
    &=C-R(\bar\theta_1\theta_1+\bar\theta_2\theta_2)
    -G(\bar\theta_1\theta_2+\bar\theta_2\theta_1)
    -D\bar\theta_1\theta_1\bar\theta_2\theta_2 \\
    &\qquad
    +(\bar\theta_1+\bar\theta_2)H
    +\bar H(\theta_1+\theta_2)
    -(\bar\theta_1\theta_1\bar\theta_2+\bar\theta_2\theta_2\bar\theta_1)\hat H
    -\bar{\hat H}(\theta_1\bar\theta_2\theta_2+\theta_1\bar\theta_1\theta_1)
  \end{aligned}
\end{equation}
and
\begin{equation}
  \mathbb M(1)
  =m+\bar\theta_1H_0+\bar H_0\theta_1-\hat m\bar\theta_1\theta_1
\end{equation}
The order parameters $C$, $R$, $G$, $D$, $m$, and $\hat m$ are ordinary numbers defined by
\begin{align}
  C=\frac{\mathbf x\cdot\mathbf x}N
  &&
  R=-i\frac{\mathbf x\cdot\hat{\mathbf x}}N
  &&
  G=\frac{\bar{\pmb\eta}\cdot\pmb\eta}N
  &&
  D=\frac{\hat{\mathbf x}\cdot\hat{\mathbf x}}N
  &&
  m=\frac{\mathbf x_0\cdot\mathbf x}N
  &&
  \hat m=-i\frac{\mathbf x_0\cdot\hat{\mathbf x}}N
\end{align}
while $\bar H$, $H$, $\bar{\hat H}$, $\hat H$, $\bar H_0$ and $H_0$ are Grassmann numbers defined by
\begin{align}
  \bar H=\frac{\bar{\pmb\eta}\cdot\mathbf x}N
  &&
  H=\frac{\pmb\eta\cdot\mathbf x}N
  &&
  \bar{\hat H}=\frac{\bar{\pmb\eta}\cdot\hat{\mathbf x}}N
  &&
  \hat H=\frac{\pmb\eta\cdot\hat{\mathbf x}}N
  &&
  \bar H_0=\frac{\bar{\pmb\eta}\cdot\mathbf x_0}N
  &&
  H_0=\frac{\pmb\eta\cdot\mathbf x_0}N
\end{align}
We can treat the integral over $\sigma_0$ immediately. It gives
\begin{equation}
  \int d\sigma_0\,e^{N\int d1\,\frac i2\sigma_0(1)(\mathbb Q(1,1)-1)}
  =2\pi\,\delta(C-1)\,\delta(G+R)\,\bar HH
\end{equation}
This therefore sets $C=1$ and $G=-R$ in the remainder of the integrand, as well
as setting everything depending on $\bar H$ and $H$ to zero. With these solutions inserted, the remaining remaining terms in the exponential give
\begin{equation}
  \begin{aligned}
    \operatorname{sdet}(\mathbb Q-\mathbb M\mathbb M^T)
    &=1+\frac{(1-m^2)D+\hat m^2-2Rm\hat m}{R^2}
    -\frac6{R^4}\bar H_0H_0\bar{\hat H}\hat H
    \\
    &\qquad+\frac2{R^3}\left[
      (mR-\hat m)(\bar{\hat H}H_0+\bar H_0\hat H)
      -(D+R^2)\bar H_0H_0
      +(1-m^2)\bar{\hat H}\hat H
    \right]
  \end{aligned}
\end{equation}
\begin{equation}
  \operatorname{sdet}f(\mathbb Q)
  =1+\frac{Df(1)}{R^2f'(1)}
  +\frac{2f(1)}{R^3f'(1)}\bar{\hat H}\hat H
\end{equation}
\begin{equation}
  \int d1\,d2\,f(\mathbb Q)^{-1}(1,2)
  =\left(f(1)+\frac{R^2f'(1)}{D}\right)^{-1}
  +2\frac{Rf'(1)}{(Df(1)+R^2f'(1))^2}\bar{\hat H}\hat H
\end{equation}

\section{Calculation of the prefactor of the average Euler characteristic}
\label{sec:prefactor}

\section{Details of the calculation of the average number of stationary points}
\label{sec:complexity.details}

\section{The quenched shattering energy in {\oldstylenums 1}\textsc{frsb} models}
\label{sec:1frsb}

\cite{Kent-DObias_2023_How}
\[
  \chi_0(q)=\frac1{\hat\omega_1}f''(q)^{-1/2}-\frac{r_d^2}{d_d}
\]

\[
  \chi(q)=\begin{cases}
    \chi_0(q) & q < q_0 \\
    1-(1-m)q_1-mq & q_0 < q < q_1 \\
    1-q & q > q_1
  \end{cases}
\]
\[
  0=\hat\omega_1r_d-\omega_1d_d
  \qquad
  \omega_1=\hat\omega_1\frac{r_d}{d_d}
\]
\[
  \log\chi
  =-\hat\omega_1 E
  +\frac12\hat\omega_1^2r_d^2/d_df'(1)
  +\frac12\int_0^1dq\,\left[
    \hat\omega_1^2f''(q)\chi(q)
    +\frac1{\chi(q)+r_d^2/d_d}
  \right]
  -\frac12\log r_d^2
\]
\[
  0=-\frac{\hat\omega_1^2f'(1)}{d_d}+\int_0^1dq\,\frac1{(r_d^2/d_d+\chi(q))^2}
\]
\[
  d_d=-\frac{1+r_d}{\int dq\,\chi(q)}r_d
\]


\bibliography{topology}

\end{document}