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

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amssymb,latexsym,graphicx}
\usepackage{newtxtext,newtxmath}
\usepackage{bbold,anyfontsize}
\usepackage[dvipsnames]{xcolor}
\usepackage[
  colorlinks=true,
  urlcolor=BlueViolet,
  citecolor=BlueViolet,
  filecolor=BlueViolet,
  linkcolor=BlueViolet
]{hyperref}

\begin{document}

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

\author{Jaron Kent-Dobias}
\email{jaron.kent-dobias@roma1.infn.it}
\affiliation{Istituto Nazionale di Fisica Nucleare, Sezione di Roma I, Rome, Italy 00184}

\begin{abstract}
  We consider the set of solutions to $M$ random polynomial equations with
  independent Gaussian coefficients on the $(N-1)$-sphere. When solutions
  exist, they form a manifold. We compute the average Euler characteristic of
  this manifold, and find different behaviors depending on the variances of the
  coefficients and $\alpha=M/N$. When $\alpha<1$, the average Euler
  characteristic is subexponential in $N$ but positive, indicating the presence
  of few connected components. When $1<\alpha<\alpha_\text{\textsc{sat}}$, it
  is exponentially large in $N$, indicating a shattering transition of the
  manifold of solutions into many components. Finally, when
  $\alpha_\text{\textsc{sat}}<\alpha$, the set of solutions vanishes. Some
  choices of variances produce $\alpha_\text{\textsc{sat}}<1$, and the shattering
  transition never takes place. We further compute the average logarithm of the
  Euler characteristic, which is representative of typical manifolds, and find
  that most of the quantitative predictions agree.
\end{abstract}

\maketitle

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,
indeed, 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.

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)=0$ for $1\leq k\leq M$. 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 $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.


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}
  \mathcal C(\mathbf x)=\frac12\sum_{k=1}^MV_k(\mathbf x)^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=\{\mathbf x\in\mathbb R^N\mid \|\mathbf x\|^2=N,V_k(\mathbf x)=0
  \;\forall\;k=1,\ldots,M\}.
\end{equation}
Because the constraints are all smooth functions, $\Omega$ is almost always a manifold without singular points. 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$ 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.

When $M$ is too large, no solutions exist and $\Omega$ becomes the empty set.
Following previous work, a replica symmetric equilibrium calculation using the
cost function \eqref{eq:cost} predicts that solutions vanish when the ratio
$\alpha=M/N$ is larger than $\alpha_\text{\textsc{sat}}=f'(1)/f(1)$. Based on the results of this paper, and the fact that this $\alpha_\text{\textsc{sat}}$ is consistent 

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=\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}
  L(\mathbf x,\pmb\omega)
  =H(\mathbf x)+\frac12\omega_0\big(\|\mathbf x\|^2-N\big)
  +\sum_{k=1}^M\omega_kV_k(\mathbf x)
\end{equation}
The integral over $\Omega$ in \eqref{eq:kac-rice} then becomes
\begin{equation} \label{eq:kac-rice.lagrange}
  \chi=\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. Details of this calculation are reserved in an appendix.

We can solve the saddle point equations in all of these
parameters save for $m=\frac1N\mathbf x_0\cdot\mathbf x$, the overlap with the
height axis. The result reduces the average Euler characteristic to
\begin{equation}
  \bar\chi\propto\int dm\,e^{N\mathcal S_\mathrm a(m)}
\end{equation}
where the annealed action $\mathcal S_a$ is given by
\begin{equation} \label{eq:ann.action}
  \begin{aligned}
    &\mathcal S_\mathrm a(m)
    =\frac12\Bigg[
      \log\left(
        \frac{\frac{f'(1)}{f(1)}(1-m^2)-1}{\alpha-1}
      \right) \\
    &\hspace{4em} -\alpha\log\left(
        \frac{\alpha}{\alpha-1}\left(
          1-\frac1{\frac{f'(1)}{f(1)}(1-m^2)}
        \right)
      \right)
    \Bigg]
  \end{aligned}
\end{equation}
and must be evaluated at a maximum with respect to $m$. This function is
plotted for a specific covariance function $f$ in Fig.~\ref{fig:action}, where
several distinct regimes can be seen.

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

  \caption{
    The annealed action $\mathcal S_\mathrm a$ of \eqref{eq:ann.action} plotted
    as a function of $m$ at several values of $\alpha$. Here, the covariance
    function is $f(q)=\frac12q^2$ and $\alpha_\text{\textsc{sat}}=2$. When
    $\alpha<1$, the action is maximized for $m^2>0$ and its value is zero. When
    $1\leq\alpha<\alpha_\text{\textsc{sat}}$, the action is maximized at
    $m=0$ and is positive. When $\alpha>\alpha_\text{\textsc{sat}}$ there is no
    maximum.
  } \label{fig:action}
\end{figure}

First, when $\alpha<1$ the action $\mathcal S_\mathrm a$ is strictly negative
and has maxima at some $m^2>0$. At these maxima, $\mathcal S_\mathrm a(m)=0$.
When $\alpha>1$, the action flips over and becomes strictly positive. In the
regime $1<\alpha<\alpha_\text{\textsc{sat}}$, there is a single maximum at
$m=0$ where the action is positive. When $\alpha\geq\alpha_\text{\textsc{sat}}$
the maximum in the action vanishes.

This results in distinctive regimes for $\overline\chi$, with an example plotted in Fig.~\ref{fig:characteristic}. If $m^*$ is the maximum of $\mathcal S_\mathrm a$, then
\begin{equation}
  \frac1N\log\overline\chi=\mathcal S_\mathrm a(m^*)
\end{equation}
When $\alpha<1$, the action evaluates to zero, and therefore $\overline\chi$ is
positive and subexponential in $N$. When $1<\alpha<\alpha_\text{\textsc{sat}}$, the action
is positive, and $\overline\chi$ is exponentially large in $N$. Finally, when
$\alpha\geq\alpha_\text{\textsc{sat}}$ the action and $\overline\chi$ are ill-defined.

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

  \caption{
    The logarithm of the average Euler characteristic $\overline\chi$ as a
    function of $\alpha$. The covariance function is $f(q)=\frac12+\frac12q^3$ and
    $\alpha_\text{\textsc{sat}}=\frac32$. The dashed line shows the average of
    $\log\chi$, the so-called quenched average, whose value differs in the
    region $1<\alpha<\alpha_\text{\textsc{sat}}$ but whose transition points
    are the same.
  } \label{fig:characteristic}
\end{figure}

We can interpret this by reasoning about topology of $\Omega$ consistent with
these results. Cartoons that depict this reasoning are shown in
Fig.~\ref{fig:cartoons}. In the regime $\alpha<1$, $\overline\chi$ is positive but not
very large. This is consistent with a solution manifold made up of few large
components, each with the topology of a hypersphere. The saddle point value
$(m^*)^2=1-\alpha/\alpha_\text{\textsc{sat}}$ for the overlap with the height axis $\mathbf x_0$ corresponds to the
latitude at which most stationary points that contribute to the Euler
characteristic are found. This means we can interpret $1-m^*$ as the typical
squared distance between a randomly selected point on the sphere and the
solution manifold.

\begin{figure}
  \includegraphics[width=0.32\columnwidth]{figs/connected.pdf}
  \hfill
  \includegraphics[width=0.32\columnwidth]{figs/shattered.pdf}
  \hfill
  \includegraphics[width=0.32\columnwidth]{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. The region of solutions is marked in black, and the critical points
    of the height function restricted to this region are marked with a point.
    For $\alpha<1$, there are few simply connected regions with most of the
    minima and maxima contributing to the Euler characteristic concentrated at
    the height $m^*$. For $\alpha\geq1$, there are many simply
    connected regions and most of their minima and maxima are concentrated at
    the equator.
  } \label{fig:cartoons}
\end{figure}

When $1<\alpha<\alpha_\text{\textsc{sat}}$, $\overline\chi$ is positive and
very large. This is consistent with a solution manifold made up of
exponentially many disconnected components, each with the topology of a
hypersphere. If this interpretation is correct, our calculation effectively
counts these components. This is a realization of a shattering transition in
the solution manifold. Here $m^*$ is zero because for any choice of height
axis, the vast majority of stationary points that contribute to the Euler
characteristic are found near the equator. Finally, for
$\alpha\geq\alpha_\text{\textsc{sat}}$, there are no longer solutions that
satisfy the constraints. The Euler characteristic is not defined for an empty
set, and in this regime the calculation yields no solution.

We have made the above discussion assuming that $\alpha_\text{\textsc{sat}}>1$.
However, this isn't necessary, and it is straightforward to produce covariance
functions $f$ where $\alpha_\text{\textsc{sat}}<1$. In this case, the picture
changes somewhat. When $\alpha_\text{\textsc{sat}}<\alpha<1$, the action
$\mathcal S_\mathrm a$ has a single maximum at $m^*=0$, where it is negative.
This corresponds to an average Euler characteristic $\overline\chi$ which is
exponentially small in $N$. Such a situation is consistent with typical
constraints leading to no solutions and a zero characteristic, but rare and
atypical configurations having some solutions.

In the regime where $\log\overline\chi$ is positive, it is possible that our
calculation yields a value which is not characteristic of typical sets of
constraints. This motivates computing $\overline{\log\chi}$, the average of
the logarithm, which should produce something characteristic of typical
samples, the so-called quenched calculation. In an appendix to this paper we
sketch the quenched calculation and report its result in the replica symmetric
approximation. This differs from the annealed calculation above only when
$f(0)>0$. The replica symmetric calculation produces the same transitions at
$\alpha=1$ and $\alpha=\alpha_\text{\textsc{sat}}$, but modifies the value
$m^*$ in the connected phase and predicts
$\frac1N\overline{\log\chi}<\frac1N\log\overline\chi$ in the shattered phase.
The fact that $\alpha_\text{\textsc{sat}}=f'(1)/f(1)$ is the same in the annealed and
replica symmetric calculations suggests that it may perhaps be exact. It is also
consistent with the full RSB calculation of \cite{Urbani_2023_A}.

We check the stability of the replica symmetric solution by calculating the
eigenvalues of the Hessian of the effective action with respect to the order
parameters. While for calculations of this kind the meaning of the sign of
these eigenvalues is difficult to understand directly, in situations where
there is a continuous \textsc{rsb} transition the sign of one of the
eigenvalues changes \cite{Kent-Dobias_2023_When}. At the $\alpha_\text{\textsc{rsb}}$ predicted in \cite{Urbani_2023_A} we see no instability of this kind, and instead only observe such an instability at $\alpha_\text{\textsc{sat}}$.

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

\begin{acknowledgements}
  JK-D is supported by a \textsc{DynSysMath} Specific Initiative of the INFN.
  The authors thank Pierfrancesco Urbani for helpful conversations on these topics.
\end{acknowledgements}

\bibliography{topology}

\paragraph{Details of the annealed calculation.}

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 H[\pmb\eta,\pmb\gamma]}
\end{align}
for real variables $\hat{\mathbf x}$ and $\hat{\pmb\omega}$, and Grassmann
variables $\bar{\pmb\eta}$, $\pmb\eta$, $\bar{\pmb\gamma}$, and $\pmb\gamma$.
With these transformations in place, there is a compact way to express $\chi$
using superspace notation. For a review of the superspace formalism for
evaluating integrals of the form \eqref{eq:kac-rice.lagrange}, see Appendices A
\& B of \cite{Kent-Dobias_2024_Conditioning}. Introducing the Grassmann indices
$\bar\theta_1$ and $\theta_1$, we define superfields
\begin{align}
  \pmb\phi(1)
  &=\mathbf x+\bar\theta_1\pmb\eta+\bar{\pmb\eta}\theta_1+\hat{\mathbf x}\bar\theta_1\theta_1 
  \label{eq:superfield.phi} \\
  \pmb\sigma(1)
  &=\pmb\omega+\bar\theta_1\pmb\gamma+\bar{\pmb\gamma}\theta_1+\hat{\pmb\omega}\bar\theta_1\theta_1
  \label{eq:superfield.sigma}
\end{align}
with which we can represent $\chi$ by
\begin{equation}
  \chi=\int d\pmb\phi\,d\pmb\sigma\,\exp\left\{
    \int d1\,L\big(\pmb\phi(1),\pmb\sigma(1)\big)
  \right\}
\end{equation}
We are now in a position to average over the distribution of constraints. Using
standard manipulations, we find the average Euler characteristic is
\begin{equation}
  \begin{aligned}
    \overline{\chi}&=\int d\pmb\phi\,d\sigma_0\,\exp\Bigg\{
      -\frac M2\log\operatorname{sdet}f\left(\frac{\phi(1)^T\phi(2)}N\right) \\
      &\qquad+\int d1\,\left[
        H\big(\phi(1)\big)+\frac12\sigma_0(1)\big(\|\phi(1)\|^2-N\big)
      \right]
    \Bigg\}
  \end{aligned}
\end{equation}

With this choice made, we can integrate over the superfields $\pmb\phi$.
Defining two order parameters $\mathbb Q(1,2)=\frac1N\phi(1)\cdot\phi(2)$ and
$\mathbb M(1)=\frac1N\phi(1)\cdot\mathbf x_0$, the result is
\begin{align}
  \overline{\chi}
  &=\int d\mathbb Q\,d\mathbb M\,d\sigma_0 \notag\\
   &\quad\times\exp\Bigg\{
    \frac N2\log\operatorname{sdet}(\mathbb Q-\mathbb M\mathbb M^T)
    -\frac M2\log\operatorname{sdet}f(\mathbb Q) \notag \\
    &\qquad+N\int d1\,\left[
      \mathbb M(1)+\frac12\sigma_0(1)\big(\mathbb Q(1,1)-1\big)
    \right]
  \Bigg\}
\end{align}
This expression is an integral of an exponential with a leading factor of $N$
over several order parameters, and is therefore in a convenient position for
evaluating at large $N$ with a saddle point. The order parameter $\mathbb Q$ is
made up of scalar products of the original integration variables in our
problem in \eqref{eq:superfield.phi}, while $\mathbb M$ contains their scalar
project with $\mathbf x_0$, and $\pmb\sigma_0$ contains $\omega_0$ and
$\hat\omega_0$. 

\paragraph{Quenched average of the Euler characteristic.}

\begin{equation}
  D=\beta R
  \qquad
  \hat\beta=-\frac{m+\sum_aR_{1a}}{\sum_aC_{1a}}
  \qquad
  \hat m=0
\end{equation}

\begin{align}
  &\mathcal S(m,C,R)
  =\frac12\log\det\big[I+\hat\beta R^{-1}(C-m^2)\big] \notag \\
  &\quad-\frac\alpha2\log\det\big[I+\hat\beta\big(R\odot f'(C)\big)^{-1}f(C)\big]
\end{align}

The quenched average of the Euler characteristic in the replica symmetric ansatz becomes for $1<\alpha<\alpha_\text{\textsc{sat}}$
\begin{align}
  \frac1N\overline{\log\chi}
  =\frac12\bigg[
    \log\left(-\frac 1{\tilde r_d}\right)
    -\alpha\log\left(
      1-\Delta f\frac{1+\tilde r_d}{f'(1)\tilde r_d}
    \right) \notag \\
    -\alpha f(0)\left(\Delta f-\frac{f'(1)\tilde r_d}{1+\tilde r_d}\right)^{-1}
  \bigg]
\end{align}
where $\Delta f=f(1)-f(0)$ and $\tilde r_d$ is given by
\begin{align}
  \tilde r_d
  =-\frac{f'(1)f(1)-\Delta f^2}{2(f'(1)-\Delta f)^2}
  \bigg(
    \alpha-2+\frac{2f'(1)f(0)}{f'(1)f(1)-\Delta f^2} \notag\\
    +\sqrt{
      \alpha^2
      -4\alpha\frac{f'(1)f(0)\Delta f\big(f'(1)-\Delta f\big)}{\big(f'(1)f(1)-\Delta f^2\big)^2}
    }
  \bigg)
\end{align}
When $\alpha\to\alpha_\text{\textsc{sat}}=f'(1)/f(1)$ from below, $\tilde r_d\to -1$, which produces $N^{-1}\overline{\log\chi}\to0$.

\end{document}