Non-Equilibrium Dynamics of Hybrid Continuous-Discrete Ground-State Sampling

## Abstract

Abstract We propose a general framework for a hybrid continuous-discrete algorithm that integrates continuous-time deterministic dynamics with Metropolis-Hastings steps to combine search dynamics with and without detailed balance. Our purpose is to study the non-equilibrium dynamics that leads to the ground state of rugged energy landscapes in this general setting. Our results show that MH-driven dynamics reach “easy” ground states faster, indicating a stronger bias in the non-equilibrium dynamics of the algorithm with reversible transition probabilities. To validate this, we construct a set of Ising problem instances with a controllable bias in the energy landscape that makes one degenerate solution more accessible than another. The constructed hybrid algorithm demonstrates significant improvements in convergence and ground-state sampling accuracy, achieving a 100x speedup on GPUs compared to simulated annealing, making it well-suited for large-scale applications.

## 1 Introduction

The fact that non-equilibrium relaxation dynamics reach different degenerate ground states at different rates, contrary to what is expected at equilibrium, is a fundamental topic in optimization and statistical physics (Bernaschi et al., ; Lucas, ; Mohseni et al., ). In machine learning, recent works suggest that sampling from more easily accessible degenerate ground states favors finding models that generalize well during the training of overparameterized neural networks (Soudry et al., ; Baity-Jesi et al., ; Feng & Tu, ; Baldassi et al., ; ). Here, we study this issue from the general perspective of non-convex optimization, focusing on the Ising Hamiltonian. In particular, we study the issue of sampling from degenerate ground states in non-convex energy landscapes, specifically the distinction between Boltzmann and chaotic sampling characterized by detailed balance (reversible) and asymmetric transition probabilities (non-reversible), respectively.

The Metropolis-Hastings (MH) algorithm is one of the most widely used techniques for sampling from a reversible Markov Chain Monte Carlo (MCMC) (Zhang et al., ; Chib & Greenberg, ; Chen et al., ). This approach proposes new states and accepts them with a probability designed to ensure that the stationary distribution of the chain is the Boltzmann distribution. Despite its simplicity and versatility, MH can struggle with getting trapped in local energy minima, especially when dealing with complex energy landscapes. Examples of these energy landscape often arise in systems exhibiting broken ergodicity (Bernaschi et al., ), such as spin glasses. In these cases, other MCMC algorithms that maintain detailed balance, like simulated annealing (Kirkpatrick et al., ) and more advanced methods (Hukushima & Nemoto, ; Surjanovic et al., ; Hukushima & Iba, ; Houdayer, ; Landau et al., ), are often favored for exploring rugged energy landscapes and finding optimal solutions. Interestingly, several algorithms do not satisfy detailed balance yet remain competitive (Leleu et al., ; ; Goto et al., ; Buesing et al., ), particularly those based on chaotic sampling (Leleu et al., ; ). Although these approaches do not sample in general from the exact Boltzmann distribution at equilibrium, it can be shown in special cases that non-reversible Markov chains mix faster than reversible ones (Kapfer & Krauth, ). However, their ability to reach degenerate ground states of rugged energy landscape during non-equilibrium dynamics is still not well understood.

In this work, we introduce a generalized algorithm that combines MCMC techniques through a MH step with chaotic search dynamics via the chaotic amplitude control (CAC) algorithm. The hyperparameters of this hybrid approach allow for a smooth transition between the two strategies. Additionally, we generate new planted instances with degenerate ground states, inspired by the Wishart planted models. A tunable parameter controls the bias towards different degenerate ground states, enabling us to investigate the non-equilibrium behavior of the hybrid MH/chaotic algorithm in locating solutions that vary in difficulty. We demonstrate that the algorithm’s performance can be optimized based on the instance bias, making it competitive with state-of-the-art methods.

## 2 Related Work

### 2.1 Discrete state MCMC optimizers

MCMC methods for optimizing energy functions over discrete spaces—such as simulated annealing  (Kirkpatrick et al., ) and parallel tempering  (Hukushima & Nemoto, )—are standard approaches in the field of optimization. Building upon these, recent schemes such as population annealing  (Hukushima & Iba, ), variational parallel tempering  (Surjanovic et al., ) and specialized techniques, including cluster updates  (Houdayer, ), Wang-Landau sampling  (Landau et al., ), and Replica Exchange Wang-Landau Sampling (Vogel et al., ), have been developed to enhance convergence and sampling efficiency. However, these methods often struggle in higher-dimensional discrete spaces due to the complexity of the energy landscape and the tendency to become trapped in local energy minima caused by broken ergodicity, which is a behavior observed in many systems with frustrated interactions. (Bernaschi et al., ).

### 2.2 Hybrid continuous-discrete samplers

For continuous space sampling, a plethora of methods have been developed based on the idea of combining continuous-time dynamics (typically via ODEs or SDEs) with discrete MCMC corrections. Some examples of this are Metropolis Adjusted Langevin Dynamics (MALA)  (Roberts & Rosenthal, ), Hamiltonian Monte Carlo (HMC)  (Neal et al., ; Xifara et al., ), as well as derivative methods such as the No-U-Turn Sampler (NUTS)  (Hoffman et al., ), Generalized HMC  (Chen et al., ), HMC with the Gaussian integral trick  (Zhang et al., ), and amortized Metropolis adjustments  (Zhang et al., ). These methods are primarily designed for continuous spaces and do not directly apply to discrete-state sampling. However, research has been conducted to extend continuous sampling techniques to discrete spaces, including methods like discontinuous HMC  (Nishimura et al., ), Gibbs with Gradient (GWG)  (Grathwohl et al., ), and the Hamming Ball Sampler  (Titsias & Yau, ). Nonetheless, these approaches are not suitable for ground-state sampling.

### 2.3 Hybrid continuous-discrete optimizers

Recently, there has been renewed interest in using continuous-space ODEs to find ground states in discrete spaces. This effort has been driven primarily by the idea of leveraging a diverse array of physical hardware, including CMOS, memristors, spintronics, superconducting circuits, and photonics (see review  (Mohseni et al., )), all designed to solve optimization problems more efficiently than current high-performance digital computers.

In particular, dynamical systems such as simulated bifurcation machines  (Goto et al., ), analog iterative machines (AIM)  (Kalinin et al., ), and chaotic amplitude control (CAC)  (Leleu et al., ; ; Reifenstein et al., ) offer an alternative approach. Rather than relying on gradient descent of an energy function in systems with symmetric interactions, such as in the seminal work by Hopfield  (Hopfield & Tank, ), it has been shown that Hamiltonian dynamics  (Goto et al., ; Neal et al., ) or chaotic dynamics driven by asymmetric interactions  (Leleu et al., ; ) can outperform standard heuristics like simulated annealing  (Kirkpatrick et al., ), parallel tempering  (Hukushima & Iba, ), and breakout local search  (Benlic & Hao, ) on specific benchmark datasets  (Leleu et al., ), challenging the prevailing belief that MCMC heuristics maintain an advantage in discrete optimization over systems whose dynamics are relaxed into the continuous domain.

Hybrid approaches that combine continuous-space ODEs with discrete MCMC probabilistic transitions remain underexplored in the context of ground-state search, as opposed to their more common application in Boltzmann sampling.

## 3 Method

### 3.1 Overview

In this section, we introduce a novel MCMC method that combines chaotic search with the MH algorithm (Hastings, ) (see Fig. 1 a). The algorithm is built on three main concepts: (1) relaxation from discrete to continuous space (or continuous embedding), (2) addition of auxiliary variables to escape from local minima (chaotic amplitude control (Leleu et al., ; ) and momentum (Kalinin et al., )), (3) probabilistic jumps in the discrete space based on the Metropolis-Hastings criterion.

In the next section, we construct a novel set of Ising problem instances inspired from the Wishart planted ensemble (Hamze et al., ) in order benchmark this algorithm’s ability to sample from optimal solutions (or ground-state energies). The constructed instances exhibit degenerate ground states that are sampled with unequal probabilities by local search algorithms (see Fig. 1 b).

Figure: Figure 1: a) Proposed hybrid approach generates deterministic chaotic trajectories that rapidly explores the energy landscape and probabilistic jumps accepted via the Metropolis-Hastings criterion. (b) Constructed planted instances with degenerate ground states with unequal sampling probabilities by heuristics. Some ground states are easier and others harder to be find by a local search heuristic.
Refer to caption: x1.png

### 3.2 Deterministic continuous time model

We define a continuous time dynamical system inspired by the CAC algorithm (Leleu et al., ; ) and generalize its formulation to draw a correspondence between its parameters and the inverse temperature $\beta$ of the Boltzmann distribution $P(\boldsymbol{\sigma})\propto e^{\beta V(\boldsymbol{\sigma})}$. The aim is to construct a Markov chain that is capable of sampling from the distribution $P(\boldsymbol{\sigma})$ defined on the discrete space $\boldsymbol{\sigma}\in\{-1,1\}^{N}$, where $N$ is the number of variables or spins.

First, discrete variables $\boldsymbol{\sigma}$ are relaxed to the continuous space $\boldsymbol{x}\in\mathbb{R}^{N}$ and we define the potential (or energy function) $V(\boldsymbol{x})$ on the real space. The following dynamical system is used to direct the search to lower energy states:

$$ $\displaystyle\gamma\frac{d^{2}\boldsymbol{u}}{dt^{2}}+\frac{d\boldsymbol{u}}{ dt}=-\alpha\boldsymbol{u}-\boldsymbol{e}\circ\nabla_{x}V,$ (1) $$

where the symbol $\circ$ denote the hadamard product; $\nabla_{x}V$, the gradient of $V$ with respect to the vector $\boldsymbol{x}$; $\alpha$, a positive parameter; and $\boldsymbol{e}$, a vector of positive auxiliary variables. The term proportional to $\gamma$ represent the momentum which has been utilized for non-convex optimization (Kalinin et al., ) and sampling (Neal et al., ; Xifara et al., ). Moreover, we set:

$$ $\displaystyle x_{i}=\phi(\tilde{\beta}u_{i}),\forall i\in\{1,\cdots,N\},$ (2) $$

and $\phi(x)=\frac{2}{1+e^{-x}}-1$ is a sigmoidal function normalized to the domain $x\in[-1,1]$.

The auxiliary variables $e_{i}$ undergo dynamic changes over time, with the goal of rectifying amplitude heterogeneity are given as follows (Leleu et al., ):

$$ $\displaystyle\frac{d\boldsymbol{e}}{dt}=-\xi(\boldsymbol{x}\circ\boldsymbol{x} -a)\circ\boldsymbol{e},$ (3) $$

with $a$ and $\xi$ being parameters: $a$ represents target amplitude of the $x$ variables and $\xi$ represents the speed of the error correction dynamics. Because this algorithm adds momentum dynamics to the original chaotic amplitude control equations (Leleu et al., ; ), we call it chaotic amplitude control with momentum (CACm).

In the Appendix section S1, the role of the auxiliary variables $\boldsymbol{e}$ is analyzed and we show that their modulation can be interpreted as a dynamic pre-conditioning of the Hessian of $V$ at the proximity of the discrete state. Because the auxiliary variables $\boldsymbol{e}$ remain positive, they do not change the inertia (number of positive, negative, and zero eigenvalues) of the Hessian of $V$ but can change its directions and eigenvalues. Thus, CAC performs some dynamic deformation of the potential function $V$ which affects its ability to sample local minima determined by the shape of $V$ when compared to other algorithms, as it will be shown in the following.

### 3.3 Probabilistic discrete time model and MCMC

In order to utilize this dynamical system algorithmically, we discretize eqs. 1-2 using a straightforward Euler approximation. Unlike Hamiltonian Monte Carlo that requires a careful selection of the integration method, an accurate integration scheme is not necessarily needed in this context, as the system can be shown to be bounded due to the influence of the sigmoid function and is dissipative. The utilization of a simple integration method and large integration time step enables the discrete-time approximation of the system to function as computationally efficient solver for combinatorial optimization. The discrete-time model is characterized as follows (after a change of variables):

$$ $\displaystyle\boldsymbol{u}(t+1)$ $\displaystyle=\boldsymbol{u}(t)-\alpha\boldsymbol{u}(t)-\boldsymbol{e}\nabla_{ x}V(\boldsymbol{x}(t))+\gamma(\boldsymbol{u}(t)-\boldsymbol{u}(t-1)),$ (4) $\displaystyle\boldsymbol{e}(t+1)$ $\displaystyle=\boldsymbol{e}(t)-\xi(\boldsymbol{x}(t)\circ\boldsymbol{x}(t)-a) \circ\boldsymbol{e}(t).$ (5) $$

To maintain the effect of the deterministic dynamics when used for optimization, we formulate an MCMC as follows. The eqs. 4 and 5 are iterated for $n$ time steps, subsequent to which we sample a discrete state $\boldsymbol{\sigma}\in\{-1,1\}^{N}$ as detailed below:

$$ $\displaystyle P(\sigma_{i}=1)=\frac{x_{i}+1}{2},\sigma_{i}=-1\text{ otherwise.}$ (6) $$

It is worth noting that all spins are updated synchronously and independently. While the deterministic dynamics over $n$ time-steps can be perceived as local updates, the generation of $\boldsymbol{\sigma}$ can leap to remote states. Note moreover that when we set $\xi=0$, $e_{i}(0)=0$, $\alpha=1$, $t=1$, and update one spin $i$ at a time, we revert to the equation of the Boltzmann machine.

### 3.4 Metropolis adjusted chaotic sampling

When $\xi\neq 0$, detailed balance is violated, causing the stationary distribution of the chain to differ from the Boltzmann distribution. To design an MCMC method that converges (at least theoretically) to the Boltzmann distribution, we introduce a Metropolis-Hastings acceptance criterion for the discrete state samples $\sigma(kt)$, where $k$ belongs to $\{1,2,...,K\}$. The inverse of the number of Metropolis-Hastings steps is denoted $f_{\text{MH}}=\frac{1}{K}$ and the total number of steps is noted $T$ with $T=Kn$.

Typically, the acceptance criterion of a new configuration $\boldsymbol{\sigma}^{2}$, when the current state is $\boldsymbol{\sigma}^{1}$ is given as follows:

$$ $\displaystyle A(\boldsymbol{\sigma}^{2}/\boldsymbol{\sigma}^{1})=\text{min}(1, \frac{P(\boldsymbol{\sigma}^{2})Q(\boldsymbol{\sigma}^{1}/\boldsymbol{\sigma}^ {2})}{P(\boldsymbol{\sigma}^{1})Q(\boldsymbol{\sigma}^{2}/\boldsymbol{\sigma}^ {1})}).$ (7) $$

Here, $Q(\boldsymbol{\sigma}^{2}/\boldsymbol{\sigma}^{1})$ represents the probability of transitioning from $\boldsymbol{\sigma}^{1}$ to $\boldsymbol{\sigma}^{2}$. The Monte Carlo step for our method is comprised of two substeps as follows: (1) A deterministic trajectory following eqs. 4 and 5 for $n$ steps starting from $\boldsymbol{u}(kn)=\vec{0}$, $\boldsymbol{e}(kn)=\vec{0}$, and , $\boldsymbol{x}(kn)=\boldsymbol{\sigma}^{1}(kn)$, (2) A probabilitic sample $\boldsymbol{\sigma}^{2}$ generated with the probability $P^{2}$ with $P^{2}=\frac{\boldsymbol{x}^{2}+1}{2}$ with $\boldsymbol{x}^{2}$ defined as $\boldsymbol{x}^{2}=\boldsymbol{x}((k+1)n)$.

It is crucial to note that $\boldsymbol{x}((k+1)n)$ is inherently reliant on $\boldsymbol{\sigma}^{1}$ due to the fact that it is generated from the path of the deterministic dynamical system which commences from $\boldsymbol{\sigma}^{1}$. Because of the synchronous updating of all spins, we can write the following definition for $Q(\sigma^{2}/\sigma^{1})$:

$$ $\displaystyle\text{log }Q(\boldsymbol{\sigma}^{2}/\boldsymbol{\sigma}^{1})= \sum_{i=1}^{N}\left[\frac{1+\boldsymbol{\sigma}^{2}_{i}}{2}\text{log }(P_{i}^{ 2})+\frac{1-\boldsymbol{\sigma}^{2}_{i}}{2}\text{log }(1-P_{i}^{2})\right].$ (8) $$

To calculate the transition probability in the reverse Monte Carlo step from $\boldsymbol{\sigma}^{2}$ to $\boldsymbol{\sigma}^{1}$, we need to apply first the deterministic sub-step (1) from $\boldsymbol{\sigma}^{2}$ to obtain a candidate end point $\tilde{\boldsymbol{x}}^{1}$, from which we can calculate the probability to generate the random sample $\boldsymbol{\sigma}^{1}$ with probability $\tilde{P}^{1}=\frac{\tilde{\boldsymbol{x}}^{1}+1}{2}$. The acceptance criterion can then be rewritten as follows:

$$ $\displaystyle\text{log }A(\boldsymbol{\sigma}^{2}/\boldsymbol{\sigma}^{1})$ $\displaystyle=\text{min}(0,\beta(V(\boldsymbol{\sigma}^{1})-V(\boldsymbol{ \sigma}^{2})-\sum_{i=1}^{N}\left[\frac{1+\boldsymbol{\sigma}^{2}_{i}}{2}\text{ log }(P_{i}^{2})+\frac{1-\boldsymbol{\sigma}^{2}_{i}}{2}\text{log}(1-P_{i}^{2} )\right]$ $\displaystyle+\sum_{i=1}^{N}\left[\frac{1+\boldsymbol{\sigma}^{1}_{i}}{2}\text {log}(\tilde{P}_{i}^{1})+\frac{1-\boldsymbol{\sigma}^{1}_{i}}{2}\text{log }(1- \tilde{P}_{i}^{1})\right]).$ (9) $$

This acceptance criterion can be computed efficiently since it only requires the logarithm of terms already computed within the deterministic equations. The discrete state $\boldsymbol{\sigma}^{2}$ is accepted with probability $A(\boldsymbol{\sigma}^{2}/\boldsymbol{\sigma}^{1})$. The Markov chain using the Metropolis-Hastings step respects detailed balance condition and converges to our target Boltzmann distribution $P(x)=e^{-\beta V(x)}$.

In summary, the pseudo code of our approach is given in section S2 of the Appendix. Note that $\boldsymbol{x}^{2}$ is solely determined by $\boldsymbol{\sigma}^{1}$. The deterministic path to $\boldsymbol{x}^{2}$ originates from $\boldsymbol{\sigma}^{1}$ with $\boldsymbol{e}=1$. Consequently, the transition probability from $\boldsymbol{\sigma}^{1}$ to $\boldsymbol{\sigma}^{2}$ depends only on $\boldsymbol{\sigma}^{1}$ (mapped to $\boldsymbol{x}^{2}$ by the deterministic flow), preserving the Markov property. The updates of all spins are done in parallel independently, resulting in a product of probabilities. This can be reduced to a sum of log probabilities for more stable numerical simulations as shown in eq. 9. We call the algorithm defined by eqs. 4-9 Metropolis Hastings Chaotic Amplitude Control with momentum (MHCACm).

### 3.5 Summary of algorithms

A key property of new algorithm proposed in this work is that can be interpreted as a generalization of many existing methods which exist as limiting cases of MHCACm. Some examples are: Simulated annealing (Kirkpatrick et al., ), Hopfield neural networks (Hopfield & Tank, ), analog iterative machines (Kalinin et al., ), chaotic amplitude control (Leleu et al., ), and chaotic amplitude control with momentum. These limiting cases are summarized in Table 1. The generalized algorithms tend to have more parameters which need to be accurately tuned in order to adequately utilize the dynamics provided by each aspect of the algorithm. When properly tuned to a given class of instances, the performance of the generalized algorithm is at least as good as all of its constituting parts. Importantly, automatic parameter tuning methods have recently been developed which allow the many parameters to be quickly and accurately tuned  (Reifenstein et al., ).

**Table 1: Limit cases of the proposed algorithm MHCACm**
| model | $f_{\text{MH}}$ | $\gamma$ | $\xi$ | model | $f_{\text{MH}}$ | $\gamma$ | $\xi$ |
| --- | --- | --- | --- | --- | --- | --- | --- |
| SA | $=\frac{1}{T}$ | $=0$ | $=0$ | CAC | $=1$ | $=0$ | $>0$ |
| HNN | $=1$ | $=0$ | $=0$ | CACm | $=1$ | $>0$ | $>0$ |
| AIM | $=1$ | $>0$ | $=0$ | MHCACm | $\frac{1}{T}<f_{\text{MH}}<1$ | $>0$ | $>0$ |

## 4 Benchmark

### 4.1 Degenerate planted Wishart instances

We construct a set of planted Ising problem instances that exhibit ground states with different tunable properties. The construction of these instances is inspired by the Wishart planted ensemble (WPE) (Hamze et al., ). WPE instances are particularly useful for benchmarking due to their tunable hardness parameter $\alpha_{\text{WPE}}$(^1^11$\alpha_{\text{WPE}}$ is set to 0.8 in this paper.). We introduce the degenerate Wishart planted ensemble (dWPE) to allow for planting of multiple degenerate ground states. Additionally, we add a tunable "bias", denoted $b$, to the energy landscape which allows us to make some ground states easier/harder to find than others for heuristic algorithms as depicted in Fig. 1 (b). Thus, we refer to these different ground states as "easy" and "hard" solutions. The details of the construction method is described in the Appendix section S5.

### 4.2 Computation time to find easy and hard solutions

In this section, we evaluate the novel algorithm for finding low energy states of non-convex combinatorial optimization problems. A common metric for evaluating the performance of Ising solvers is the "time to solution" (TTS) which measures the number of steps needed to have 99$\%$ probability of finding any ground state. In case there are $D$ degenerate ground states (i.e., solutions of same energy), there is equal probability of sampling them from the Boltzmann distribution. However, algorithms generally reach a ground state before equilibration to the Boltzmann distribution and, consequently, each degenerate solution can have different probability of being sampled at finite sampling time. We define TTS, $\text{TTS}_{\text{easy}}$, and $\text{TTS}_{\text{hard}}$ as the time to reach any, the easier, and harder to find degenerate solution (of same Ising energy) with 99$\%$ probability as follows:

$$ $\displaystyle\text{TTS}=T\frac{\log(1-0.99)}{\log\left(1-\sum_{d}P_{d}\right)} ,\text{TTS}_{\text{easy}}=T\frac{\log(1-0.99)}{\log\left(1-\text{min}_{d}P_{d} \right)},\text{TTS}_{\text{hard}}=T\frac{\log(1-0.99)}{\log\left(1-\text{max}_ {d}P_{d}\right)},$ (10) $$

where $P_{d}$ is the probability to sample the ground state $d$. The algorithms used for comparison are all tuned to near optimal parameters by minimizing the TTS. We use a state-of-the-art autotuning method called dynamic anisotropic smoothing (Reifenstein et al., ) (DAS) which is particularly well fitted to optimizing the parameters of heuristics for combinatorial optimization  (Reifenstein et al., ). The use of a automatic parameter tuning also guarantees a fair comparison between algorithms (see Appendix S4). Moreover, all algorithms are compared using the number of “steps" to solution, where the algorithmic complexity of a single step is dominated by the matrix-vector multiplication of size $N\times N$(^2^22For SA, a step represents a full sweep of $N$ spin updates which also has complexity $O(N^{2})$. The same for PT, but divided by the number of replica.).

### 4.3 Global optimization to any optimal solution

For the task of finding any ground state, Fig. 2 (a) shows that MHCACm performs almost as well as CACm and better than AIM, CAC, and SA(^3^33Results for SA and PT are collected using dwave-neal (dwavesystems, ) and pySA (Mandrà & Munoz-Bauza, )). In the case of sampling any ground state of biased ($b=12$ ) dWPE problem instances (see section 4.1), MHCACm requires significantly less steps than other algorithm to find an optimal solution (see Fig. 2 (b)).

Figure: Figure 2: The time to solution of (unbiased) Wishart planted instances and biased degenerate planted instance with $b=12$ are shown in (a) and (b), respectively. $f_{\text{MH}}=0.1$. All other parameters tuned with DAS (Reifenstein et al., ) (see Appendix S4). Colored regions represent the 95% confidence interval calculated by generating 15 instances per data point.
Refer to caption: x2.png

### 4.4 Sampling in the case of biased degenerate optimal solutions

Next, we examine why MHCACm requires fewer steps to find any solution when there is a larger bias in the dWPE instances. Fig. 3 (a) and (b) show that the $\text{TTS}_{\text{easy}}$ for easy ground states decreases as a function of the bias $b$ of the generated Wishart instances. Conversely, the time to sampled the hard solution $\text{TTS}_{\text{hard}}$ increases. This is unlike the case of CACm without the MH step, for which the $\text{TTS}_{\text{easy}}$ and $\text{TTS}_{\text{hard}}$ remain approximately equal for all bias $b>0$. Thus, MHCACm is better at finding the easier solutions, which explains the smaller overall TTS (to any ground state) observed in Fig. 2. Table 2 indicates that MHCACm is the fastest algorithm to sample from easier ground states, whereas CACm is the better algorithm to reach the harder to find ones.
This difference in complexity between finding the easy and hard solutions is also present in other MCMC algorithms such as SA and PT as can be seen in Tab. 2. The energy landscape is structured so that there are many more low energy states in the vicinity of the easy solution. We also show the results of the Boltzmann sampling algorithm Gibbs with gradient (Grathwohl et al., ) (GWG), although its purpose is not combinatorial optimization. An algorithm that attempts to sample from a Boltzmann distribution will thus have a large bias towards this part of the solution space at finite temperature, because of the many low energy states that should be sampled with high probability. Although in the infinite time limit we would expect equal probability of finding all ground states, in practice the extremely slow relaxation time of these algorithms makes it very improbable for the "hard" ground state to be found. However, our hybrid algorithm is able to find both ground states with nonzero probability. We believe a reason for this is that the purely analog methods such as CAC, AIM and CACm do not sample from a Boltzmann distribution and thus can partially avoid these large basins of attraction. We see this in Tab. 2 in which the analog dynamics are able to find both ground states with reasonable success as well.

Figure: Figure 3: Time to sample the easy solution $\text{TTS}_{\text{easy}}$ (a) and hard solution $\text{TTS}_{\text{hard}}$ (b) vs. the bias $b$ of degenerate Wishart planted instances. (c) $\text{TTS}_{\text{easy}}$ vs. runtime $T$ for degenerate Wishart planted instances. The same optimal parameters as in Fig. 2 are used. (a,b) $N=100$. (c) $N=100$, $b=10$.
Refer to caption: x3.png

**Table 2: Time to solution (in steps, i.e., matrix vector multiplication) using degenerate Wishart planted instances with bias $b=12$. We show both TTS, $\text{TTS}_{\text{easy}}$ and $\text{TTS}_{\text{hard}}$ (as defined in the text) obtained for a variety of algorithms and two problem sizes. Since SA and GWG were not able to sample the hard ground state, $\text{TTS}_{\text{hard}}$ is not defined for it and $\text{TTS}_{\text{hard}}$ = TTS. Bold numbers denote the best algorithm.**
|  | $10^{-3}$ TTS | $10^{-3}$ $\text{TTS}_{\text{easy}}$ | $10^{-3}$<br>$\text{TTS}_{\text{hard}}$ |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| model/N | $100$ | $140$ | $100$ | $140$ | $100$ | $140$ |
| SA | 27.56 | 1917.67 | 27.56 | 1917.67 | N/A | N/A |
| PT | 2437.92 | N/A | 287.71 | N/A | 4605.06 | N/A |
| GWG | 6369.48 | 15779.94 | 6369.48 | 15779.94 | N/A | N/A |
| AIM | 22.31 | 906.96 | 37.06 | 958.13 | 390.72 | 16884.93 |
| CAC | 238.81 | 5833.26 | 586.26 | 46960.90 | 427.34 | 6634.85 |
| CACm | 9.41 | 22.91 | 18.54 | 37.84 | 25.08 | 76.88 |
| MHCACm | 1.40 | 2.44 | 1.47 | 2.44 | 318.16 | 248.63 |

### 4.5 Results on GSET

The GSET  (Ye, ) is a set of Max-Cut problem instance which are commonly used to benchmark Ising and Max-Cut solvers  (Goto et al., ; Leleu et al., ; Benlic & Hao, ). To test the ability of MHCACm to find optimal solutions on larger optimization problems, we briefly present some results on these instances in comparison with a state of the art algorithm known as dSBM (Goto et al., ). Table 3 shows that MHCACm reaches best known solutions faster than dSBM (Goto et al., ) on the first few instances of the GSET problem set (Ye, ).

**Table 3: TTS (in steps or MVM) for solving GSET (Ye, ) instances of MHCACm and dSBM (Goto et al., ).**
| instance | MHCACm | dSBM |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| $N$ (id) | $p_{0}$ | T | TTS | $p_{0}$ | T | TTS |
| 800 (1) | 0.612000 | 3000 | 14592 | 0.026777 | 2000 | 339332 |
| 800 (2) | 0.021000 | 3000 | 650949 | 0.010660 | 6000 | 2578124 |
| 800 (3) | 0.533500 | 3000 | 18118 | 0.033920 | 3000 | 400343 |
| 800 (4) | 0.140500 | 3000 | 91249 | 0.025144 | 2000 | 361673 |
| 800 (5) | 0.044000 | 3000 | 307029 | 0.022099 | 3000 | 618221 |
| 800 (6) | 0.070500 | 3000 | 188972 | 0.023856 | 1000 | 190728 |
| 800 (7) | 0.300000 | 3000 | 38734 | 0.022552 | 1000 | 201889 |
| 800 (8) | 0.166500 | 3000 | 75858 | 0.019060 | 1500 | 358947 |

### 4.6 Parallelization on GPU

The algorithm is particularly well-suited for deployment on highly parallel hardware, such as GPUs, due to its capacity for parallel spin updates. In particular, the computational bottleneck of MHCACm is the matrix vector multiplication (MVM), a basic linear algebra operation that can be done very efficiently on modern GPUs using frameworks such as PyTorch. In table 4 we compare the wall-clock TTS of our PyTorch implementation of MHCACm finding ground states of a degenerate WPE instance of size $N=100$. We also show the wall-clock time of a standard implementation of SA  (Tiosue, ) running an CPU. MHCACm on GPU exhibits a 100x speed up in real time against simulated annealing on CPU for sampling easy optimal solutions. Furthermore, MHCACm is compatible with an implementation on special-purpose analog hardware designed for linear algebra acceleration, such as a memresistor crossbar (Hu et al., ) or on-chip photonic mesh arrays (Bogaerts et al., ).

**Table 4: Wall clock times of algorithms on different computing platforms. TTS is time to find any ground state of a degenerate WPE instance with $\alpha_{\text{WPE}}=0.8,N=100$ and bias $b$ = 12. The GPU used is a Nvidia V100 and the CPU wall time is calculated on a 2024 MacBook Pro (Apple M3 Chip).**
|  | TTS (s) ($N=100$) |  |
| --- | --- | --- |
|  | unbiased WPE | biased WPE ($b=12$) |
| SA (CPU) | $5.794\times 10^{-7}$ | $0.0217$ |
| MHCACm (CPU) | $1.739\times 10^{-6}$ | $0.00243$ |
| MHCACm (GPU) | $\mathbf{9.798\times 10^{-8}}$ | $\mathbf{1.373\times 10^{-4}}$ |

## 5 Comparison to other approaches

Methods such as HMC (Neal et al., ; Xifara et al., ) introduce auxiliary degrees of freedom (momentum) to the gradient dynamics, providing the system additional directions to move away from local minima of the error function. However, these methods do not force the dynamics to be close to the discrete state, rendering them ill-suited for discrete state sampling. CAC suggests a different mechanism for the auxiliary variables. It adjusts the amplitude of the continuous variables, ensuring proximity to the discrete (binary) state (Leleu et al., ; ). Conventionally, this can be interpreted as a dual-primal Lagrangian problem in which we want to minimize the error function (in the continuous embedding) while also satisfying a constraint (being close to the discrete state). Our method differs from the primal dual (Vadlamani et al., ) approach. With CAC, the error function is achieved by multiplying the gradient with auxiliary variables, akin to dynamically pre-conditioning the Hessian of the cost function.

Unlike HMC however, CAC does not conserve (even approximately) energy. In HMC (Neal et al., ; Xifara et al., ), the MH step compensates for the discretization errors brought about by the leapfrog integration. However, the reason for adding an MH step to CAC is different. State-of-the-art algorithms that exploit relaxation to an continuous state form discrete combinatorial optimization such as memcomputing (Sheldon et al., ) chaotic amplitude control (Leleu et al., ), and simulated bifurcation machine (Goto et al., ) do not sample fairly from the Boltzmann distribution. Our proposed methodology bridges this gap, integrating continuous relaxation algorithms with fair sampling of a discrete distribution via the construction of a hybrid MCMC.

## 6 Conclusion

The MHCACm algorithm effectively combines chaotic dynamics with the Metropolis-Hastings method to enhance ground-state sampling efficiency in discrete spaces, ensuring fair sampling at equilibrium while biasing towards more reachable ground states during non-equilibrium dynamics. This approach significantly improves the speed and accuracy of finding optimal solutions for combinatorial optimization tasks, demonstrating superior performance and suitability for large-scale parallel processing on GPUs.

The theoretical understanding of the sampling capabilities of MHCACm are out of the scope of this paper. However, numerical evidence seems to indicate that the algorithm not only achieves fair sampling from the Boltzmann distribution after sufficient relaxation time (see Apendix Fig. S2) but can exhibit faster relaxation times on NP-hard Ising problem instances with complicated energy landscape (see Fig. 3 and Tab. 1) when compared to traditional methods such as SA. Lastly, the connection with learning in over-parameterized neural networks is to be investigated. Finding the most reachable and common solution to a large scale optimization problem is particularly useful for training neural networks in machine learning, as it typically corresponds to a trained network with better generalization properties (Soudry et al., ; Baity-Jesi et al., ; Feng & Tu, ; Baldassi et al., ; ).

## - Appendix - Non-Equilibrium Dynamics of Hybrid Continuous-Discrete Ground-State Sampling

## S1 Analysis of the deterministic path

In the following section, we analyze the deterministic dynamical system defined as chaotic amplitude control with momentum (CACm) in the main manuscript. For simplicity, we do not take into account the effect of momentum and set $\gamma=0$.

### S1.1 Dynamic pre-conditioning of the Hessian of V 𝑉 V italic_V

The equation (1) of the main manuscript is a straightforward gradient descent of the potential $V$. However, the gradient is modulated by the auxiliary variables. In this section, we give more context about the motivation for using the auxiliary variables $\boldsymbol{e}$ and explain that their modulation can be interpreted as a dynamic pre-conditioning of the Hessian of $V$ at the proximity of the discrete state $\boldsymbol{\sigma}$. First, we show in Appendix section S1.2 that the variables $\boldsymbol{u}$ in the system of eqs. (1) and (2) can be eliminated by a change of variables. Then, the Jacobian matrix $J_{xx}$ in the space $\boldsymbol{x}$ defined as $J_{xx}=\{\frac{\partial\dot{x}_{i}}{\partial x_{j}}\}_{ij}$ with $\dot{x}_{i}=\frac{d\boldsymbol{x}}{dt}$ can be expressed as follows at proximity of the fixed points of the system (see Appendix section S1.3):

$$ $\displaystyle J_{xx}=D[-\alpha f(\sqrt{a}\boldsymbol{\sigma})]-D[\boldsymbol{e }\circ g(\sqrt{a}\boldsymbol{\sigma})]-D[\boldsymbol{e}h(\sqrt{a}\boldsymbol{ \sigma})]H(\sqrt{a}\boldsymbol{\sigma})$ (S1) $$

where $D[\boldsymbol{x}]$ is a diagonal matrix with elements given by the components of the vector $\boldsymbol{x}$ and $H(\sqrt{a}\boldsymbol{\sigma})$ is the Hessian of $V$ defined as the Jacobian of $\nabla_{x}V$ and calculated at the point $\sqrt{a}\boldsymbol{\sigma}$. The functions $f$ and $g$ and defined in Appendix section S1.2. Equation S1 shows that the auxiliary variables $\boldsymbol{e}$ modulates the Hessian $H$ on the slow time scale such that $\tilde{H}(t,\tau)=D[e(\tau)h(x_{i})]H$ with $h(x)=\frac{1}{\frac{d\psi}{dx}}>0$ for $x\in]-1,1[$ and $\psi(x)=\frac{\phi^{-1}(x)}{\tilde{\beta}}$ (see lemma 1 in Appendix section S1.3).

Given that the fixed points of eqs. (1) and (2) are given as $e_{i}=-\alpha\sqrt{a}\frac{1}{\sigma_{i}\frac{\partial V}{\partial x_{i}}(
\sqrt{a}\sigma_{i})}$, we can interpret the role of $\boldsymbol{e}$ as a normalization of each row of the Hessian by the inverse of the gradient in the corresponding direction (see theorem 1 in Appendix section S1.3). Because the auxiliary variables $\boldsymbol{e}$ remain positive, they do not change the inertia (number of positive, negative, and zero eigenvalues) of the Hessian $H$ (see lemma 2 in Appendix section S1.4) but can change its directions and eigenvalues. Moreover, we show in the theorem 2 of in Appendix section S1.5 that this dynamical system can exhibit period cycles and chaos when the target amplitude $a$ is chosen to be sufficiently small. In summary, the role of the auxiliary variables is to act as a pre-conditioning of the Hessian and introduce entropy in the search dynamics due to its chaotic dynamics. Similar dynamics have been shown to be competitive with state-of-the-art heuristics to solve Ising problems (Leleu et al., ; ).

### S1.2 Change of variable

Consider $x_{i}=\phi(2\beta u_{i})\implies u_{i}=\psi_{\beta}(x_{i})$ with $\psi_{\beta}(x_{i})=\frac{\phi^{-1}(x_{i})}{2\beta}$.

In the case $\phi(x)=\frac{2}{1+e^{-x}}-1$, we have $\phi^{-1}(x)=-ln(\frac{2}{x+1}-1)$, $\psi_{\beta}^{\prime}(x)=-\frac{1}{\beta(x+1)(x-1)}$ with $\psi_{\beta}^{\prime}(x)=\frac{\partial\psi_{\beta}(x)}{\partial x}$. Moreover, $\psi_{\beta}^{\prime}(x)$ is positive on $]-1,1[$ and $\frac{1}{\psi_{\beta}^{\prime}(x)}\rightarrow 0$ for $|x|\rightarrow 0$.

Because $\frac{du_{i}}{dt}=\psi_{\beta}^{\prime}(x_{i})\frac{dx_{i}}{dt}$ with $\psi_{\beta}^{\prime}=\frac{1}{\beta}\frac{\partial\phi^{-1}(x_{i})}{\partial x
_{i}}$, the equations of motion can be rewritten by eliminating $u_{i}$ ($\forall i\in\{1,\cdots,N\}$) as follows:

$$ $\displaystyle\frac{dx_{i}}{dt}$ $\displaystyle=-\alpha\frac{\psi_{\beta}(x_{i})}{\psi_{\beta}^{\prime}(x_{i})}- \frac{e_{i}}{\psi_{\beta}^{\prime}(x_{i})}\nabla_{x}V,$ (S2) $\displaystyle\frac{de_{i}}{dt}$ $\displaystyle=-\xi(x_{i}^{2}-a)e_{i},$ (S3) $$

### S1.3 Jacobian matrix at the fixed points

The fixed points $\frac{dx_{i}}{dt}=0$ and $\frac{de_{i}}{dt}=0$ of the dynamical system described by eqs. S2 and S3 are given as follows:

$$ $\displaystyle x_{i}$ $\displaystyle=\sigma_{i}\sqrt{a},$ (S4) $\displaystyle e_{i}$ $\displaystyle=-\frac{\alpha\psi_{\beta}(\sigma_{i}\sqrt{a})}{\frac{\partial V} {\partial x_{i}}(\sigma\sqrt{a})}.$ (S5) $$

Moreover, the Jacobian matrix $J$ at the fixed points is defined as:

$$ $J=\left[\begin{array}[]{cc}J_{xx}&J_{xe}\\ J_{ex}&J_{ee}\end{array}\right]$ (S6) $$

with

$$ $\displaystyle J_{xx}$ $\displaystyle=D[-\alpha f(\sqrt{a}\boldsymbol{\sigma})]-D[\boldsymbol{e}\circ g (\sqrt{a}\boldsymbol{\sigma})]-D[\boldsymbol{e}\circ h(\sqrt{a}\boldsymbol{ \sigma})]H(\sqrt{a}\boldsymbol{\sigma}),$ (S7) $\displaystyle J_{xe}$ $\displaystyle=-D[\{\frac{\frac{\partial V}{\partial x_{i}}(\sqrt{a}\sigma_{i}) }{\psi_{\beta}^{\prime}(\sqrt{a}\sigma_{i})}\}_{i}],$ (S8) $\displaystyle J_{ex}$ $\displaystyle=-2\xi\sqrt{a}D[\boldsymbol{\sigma}\circ\boldsymbol{e}],$ (S9) $\displaystyle J_{ee}$ $\displaystyle=0,$ (S10) $$

with $H(x)=J[\nabla_{x}V]$ is the Hessian of $V$ in $x$ without its diagonal elements. The notation $D[\boldsymbol{x}]$ denote the diagonal matrix with diagonal elements given by the components of the vector $\boldsymbol{x}$. The functions $f$, $g$, and $h$ of eq. S7 are defined as follows:

$$ $\displaystyle f(x_{i})$ $\displaystyle=\frac{\partial\frac{\psi_{\beta}}{\psi_{\beta}^{\prime}(x_{i})}} {\partial x_{i}},$ (S11) $\displaystyle g(x_{i})$ $\displaystyle=\frac{\partial\frac{1}{\psi_{\beta}^{\prime}(x_{i})}\frac{ \partial V}{\partial x_{i}}}{\partial x_{i}},$ (S12) $\displaystyle h(x_{i})$ $\displaystyle=\frac{1}{\psi_{\beta}^{\prime}(x_{i})}>0\text{ for }x_{i}\in]-1, 1[.$ (S13) $$

In the case of $\phi(x)=\frac{2}{1+e^{-x}}-1$, we have:

$$ $\displaystyle f(x_{i})$ $\displaystyle=2+x_{i}\text{ln}(-\frac{-1+x_{i}}{1+x_{i}}),$ (S14) $\displaystyle g(x_{i})$ $\displaystyle=-\beta x_{i}\frac{\partial V}{\partial x_{i}}-\beta(x_{i}+1)(x_{ i}-1)\frac{\partial^{2}V}{\partial x_{i}^{2}},$ (S15) $\displaystyle h(x_{i})$ $\displaystyle=-\beta(x_{i}+1)(x_{i}-1).$ (S16) $$

Note that $f$ is an even function with $f(x)=f(-x)$.

In the case $\xi\rightarrow 0$, the dynamics of $e_{i}$ variables becomes much slower than that of $x_{i}$. We note $\tau$ the slower time scale of $e_{i}$ with $\tau=t\xi$. The lemma 1 is obtained using eq. S7.

Lemma 1: When $\xi\rightarrow 0$, the Hessian $H$ of $V$ is modulated on a slow time scale $\tau$ with $\tilde{H}=D[\boldsymbol{e}\circ g(\sqrt{a}\boldsymbol{\sigma})]H$ as proximity of the fixed points $\boldsymbol{x}=\sqrt{a}\boldsymbol{\sigma}$. We call $\tilde{H}$ the effective Hessian matrix.

### S1.4 Effective Hessian

Stability of the fixed points is determined by the eigenvalues of the Jacobian matrix $J$. The blocks $J_{ex}$ and $J_{xe}$ are such that (see also (Leleu et al., )):

$$ $\displaystyle J_{ex}J_{xe}=bI,$ (S17) $$

with $I$ the identity matrix and $b=-2\xi\alpha\sqrt{a}\frac{\psi_{\beta}(\sqrt{a})}{\psi_{\beta}^{\prime}(\sqrt
{a})}$ and $J_{ee}=0$. Consequently, the determinant $P(\lambda)=|J-\lambda I|$ characterizing the eigenvalues of $J$ has the following property:

$$ $\displaystyle P(\lambda)=\left|\begin{array}[]{cc}J_{xx}-\lambda I&J_{xe}\\ J_{ex}&-\lambda I\end{array}\right|=|J_{xx}-\lambda I|\lambda-J_{xe}J_{ex}=|J_ {xx}-\lambda I|\lambda-bI.$ (S20) $$

If we note $\mu_{i}$ the eigenvalues of $J_{xx}$, the characteristic polynomial can be rewritten:

$$ $\displaystyle P(\lambda)=-(\mu-\lambda)\lambda-bI=\lambda^{2}-\mu\lambda-b.$ (S21) $$

The eigenvalues $\lambda_{j}$ can come in pairs $\lambda_{j}^{+}$ and $\lambda_{j}^{-}$ depending on the sign of $\Delta_{j}=\mu_{j}^{2}+4b$ and are expressed as follows by solving $P(\lambda)=0$ ($\forall j\in\{1,\cdots,N\}$):

$$ $\displaystyle\lambda_{j}^{\pm}=(\mu_{j}\pm\sqrt{\Delta_{j}})\frac{1}{2}\text{ if }\Delta_{j}>0,$ (S22) $\displaystyle\lambda_{j}^{\pm}=(\mu_{j}\pm i\sqrt{\Delta_{j}})\frac{1}{2}\text { if }\Delta_{j}<0.$ (S23) $$

We define the pre-conditioned Hessian matrix $\tilde{\tilde{H}}$ as $\tilde{\tilde{H}}(\sqrt{a}\boldsymbol{\sigma})=D[\frac{\sigma_{i}}{\frac{
\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})}]H(\sqrt{a}\sigma_{i})$ such that :

$$ $\displaystyle D[e_{i}h(\sqrt{a}\sigma_{i})]H(\sqrt{a}\boldsymbol{\sigma})$ $\displaystyle=D[-\alpha\frac{1}{\frac{\partial V}{\partial x_{i}}(\sigma_{i} \sqrt{a})}\frac{\psi_{\beta}(\sigma_{i}\sqrt{a})}{\psi_{\beta}^{\prime}(\sigma _{i}\sqrt{a})}]H(\sqrt{a}\boldsymbol{\sigma}),$ (S24) $\displaystyle=-\alpha\frac{\psi_{\beta}(\sqrt{a})}{\psi_{\beta}^{\prime}(\sqrt {a})}D[\frac{\sigma_{i}}{\frac{\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a}) }]H(\sqrt{a}\boldsymbol{\sigma}),$ (S25) $\displaystyle=-\alpha\frac{\psi_{\beta}(\sqrt{a})}{\psi_{\beta}^{\prime}(\sqrt {a})}\tilde{\tilde{H}}(\sqrt{a}\boldsymbol{\sigma}),$ (S26) $$

The last step uses the fact that $\frac{\psi_{\beta}(x)}{\psi_{\beta}^{\prime}(x)}$ is an odd function.

Note that the factors $\frac{\sigma_{i}}{\frac{\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})}$, $\forall i$, multiplying the Hessian $H$ are positive at local minima satisfying $\frac{\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})\sigma_{i}>0$, $\forall i$. In this case, the eigenvalues of $D[\frac{\sigma_{i}}{\frac{\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})}]H(
\sqrt{a}\sigma_{i})$ are the same as $D[\frac{\sigma_{i}}{\frac{\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})}]^{1
/2}H(\sqrt{a}\sigma_{i})D[\frac{\sigma_{i}}{\frac{\partial V}{\partial x_{i}}(
\sigma_{i}\sqrt{a})}]^{1/2}$. By Sylvester’s law of matrix inertia, this implies that the eigenvalues of $\tilde{\tilde{H}}$ and $H$ have the same signs. In matrix theory, the term ”inertia" refers to the number of positive, negative, and zero eigenvalues of a real symmetric matrix.

Lemma 2: The signs of the eigenvalues of $\tilde{\tilde{H}}$ and $H$ are the same when $\sigma_{i}\frac{\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})>0$, $\forall i$.

Theorem 1: The effect of error variables at proximity of fixed points is to replace the Hessian $H$ of $V$ by an effective Hessian $\tilde{\tilde{H}}$ with $\tilde{\tilde{H}}(\sqrt{a}\boldsymbol{\sigma})=D[\frac{\sigma_{i}}{\frac{
\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})}]H(\sqrt{a}\sigma_{i})$. The inertia of the Hessian and effective Hessian are the same for local minima of $V$ verifying the condition $\sigma_{i}\frac{\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})>0$, $\forall i$.

### S1.5 Eigenvalues of H 𝐻 H italic_H and stability of the fixed points

The eigenvalues of the Jacobian $J_{xx}$, noted $\mu_{j}$, can be expressed using the eigen spectrum of the pre-conditioned Hessian $\tilde{\tilde{H}}$, noted $\gamma_{j}$, as follows using eq. S7:

$$ $\displaystyle\mu_{j}=-\alpha f(\sqrt{a}\sigma_{i})-e_{i}g(\sqrt{a}\sigma_{i})+ \alpha\frac{\psi_{\beta}(\sqrt{a})}{\psi_{\beta}^{\prime}(\sqrt{a})}\gamma_{j}.$ (S27) $$

The fixed points become stable when the real parts of eigenvalues are all negative, which is verified when the eigenvalue with maximal real part becomes equal to 0 given by the following condition using eqs. S22 and S23:

$$ $\displaystyle\text{max}_{j}(\text{Re}[\mu_{j}\pm\sqrt{\mu_{j}^{2}+4b}])=0\text { if }\mu_{j}^{2}-4b>0,$ (S28) $\displaystyle\text{max}_{j}(\text{Re}[\mu_{j}\pm i\sqrt{-\mu_{j}^{2}-4b}])=0 \text{ if }\mu_{j}^{2}-4b<0.$ (S29) $$

Noting $\tilde{\gamma}$ the eigenvalue with maximal part $\text{Re}[\gamma_{j}]$, eqs. S28 and S28 imply that the fixed point $\sigma_{i}\sqrt{a}$ becomes stable under the condition (see eq. S23):

$$ $\displaystyle\tilde{\mu}=0,$ (S30) $\displaystyle\text{i.e., }-\alpha f(\sqrt{a}\sigma_{i})-e_{i}g(\sqrt{a}\sigma_ {i})-\tilde{\gamma}=0.$ (S31) $$

Assuming $\frac{\partial^{2}V}{\partial x_{i}^{2}}=0$, we have (see eqs. S3 and S12):

$$ $\displaystyle e_{i}g(x_{i})$ $\displaystyle=\frac{\psi_{\beta}^{\prime\prime}(x_{i})}{(\psi_{\beta}^{\prime} (x_{i}))^{2}}\frac{\partial V}{\partial x_{i}}\frac{\alpha\psi_{\beta}(x_{i})} {\frac{\partial V}{\partial x_{i}}},$ (S32) $\displaystyle=\alpha\frac{\psi_{\beta}^{\prime\prime}(x_{i})}{(\psi_{\beta}^{ \prime}(x_{i}))^{2}}\psi_{\beta}(x_{i}).$ (S33) $$

Moreover, we have (see eq. S11):

$$ $\displaystyle f(x_{i})$ $\displaystyle=\frac{\partial\frac{\psi_{\beta}}{\psi_{\beta}^{\prime}(x_{i})}} {\partial x_{i}}$ (S34) $\displaystyle=2-\frac{\psi_{\beta}\psi_{\beta}^{\prime\prime}}{(\psi_{\beta}^{ \prime}(x_{i}))^{2}}.$ (S35) $$

Combining eqs. S33 and S35 , eq. S31 can be written as follows:

$$ $\displaystyle\tilde{\gamma}(\sqrt{a}\boldsymbol{\sigma})=2\frac{\psi_{\beta}^{ \prime}(\sqrt{a})}{\psi_{\beta}(\sqrt{a})}.$ (S36) $$

Theorem 2: For small values of the target amplitude $a$ such that $\tilde{\gamma}(\sqrt{a}\boldsymbol{\sigma})>2\frac{\psi_{\beta}^{\prime}(\sqrt
{a})}{\psi_{\beta}(\sqrt{a})}$ where $\tilde{\gamma}(\sqrt{a}\boldsymbol{\sigma})$ is the eigenvalue of $\tilde{\tilde{H}}(\sqrt{a}\boldsymbol{\sigma})$ with $\tilde{\tilde{H}}(\sqrt{a}\boldsymbol{\sigma})=D[\frac{\sigma_{i}}{\frac{
\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})}]H(\sqrt{a}\sigma_{i})$, the system described by eqs. S2 and S3 does not have any stable fixed points assuming $\frac{\partial^{2}V}{\partial x_{i}^{2}}=0$.

In general settings where the system is bounded, Theorem 2 implies that the dynamical system of eqs. S2 and S3 settles to either a periodic (quasi periodic) or chaotic attractor for sufficiently large runtime T.

In the case of the Ising Hamiltonian $V(\boldsymbol{x})=-\frac{1}{2}\sum_{ij}\omega_{ij}x_{j}x_{i}$, the calculations are simplified because $\frac{\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})=\sqrt{a}\frac{\partial V
}{\partial x_{i}}(\sigma_{i})$ and $V(\boldsymbol{\sigma}\sqrt{a})=aV(\boldsymbol{\sigma})$ and eq. S36 can be written:

$$ $\displaystyle\tilde{\gamma}(\boldsymbol{\sigma})=2\sqrt{a}\frac{\psi_{\beta}^{ \prime}(\sqrt{a})}{\psi_{\beta}(\sqrt{a})},$ (S37) $$

where $\tilde{\gamma}$ is independent of $\sqrt{a}$.

In Fig. S1, the stability of the dynamics described in eqs. S2 and S3 when $V$ is Ising Hamiltonian with $N=8$ is shown by the bifurcation diagram in the space $\{F(\sqrt{a}),\sqrt{a}\}$. The $+$ and $\times$ markers correspond to the position real part of the eigenvalue of $\tilde{\gamma}$ of $D[\frac{\sigma_{i}}{\frac{\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})}]H(
\sqrt{a}\sigma_{i})$ for the x-axis and the value of $a=0.6$ and $a=0.8$, respectively. For $a=0.8$, the point $\{\tilde{\gamma},\sqrt{a}\}$ is above the set defined by $\{\sqrt{a},F(\sqrt{a})\}$ with $a>0$ (shown in red). In this regime, the system settles to a stable fixed point as shown in Fig. S1 (b) and (d). For $a=0.6$, the point $\{\tilde{\gamma},\sqrt{a}\}$ is below the red curve. In this regime, the system does not possess any fixed points and the limit behavior is a periodic cycle.

Figure: Figure S1: (a) Bifurcation diagram of the system of eqs. S2 and S3 in the case of the $V(\boldsymbol{x})=-\frac{1}{2}\sum_{ij}\omega_{ij}\sigma_{i}\sigma_{j}$ with $N=8$ and $\omega_{ij}$ is randomly chosen to be equal to $1$ or $-1$ for $j>i$. $\omega_{ij}=\omega_{ji}$ and $\omega_{ii}=0$. The red curve corresponds to the set $\{\sqrt{a},F(\sqrt{a})\}$ with $a>0$ which denote the condition under which the real part of the eigenvalues $\mu$ of the Jacobian matrix $J_{xx}$ is zero. The $+$ and $\times$ markers correspond to the maximum real part $\tilde{\gamma}$ of the eigenvalues of $D[\frac{\sigma_{i}}{\frac{\partial V}{\partial x_{i}}(\sigma_{i}\sqrt{a})}]H( \sqrt{a}\sigma_{i})$ for the x-axis and the value of $a=0.6$ and $a=0.8$, respectively. (b) and (c) Time series of soft spins $x_{i}$ vs. time in the case $a=0.8$ and $a=0.6$, respectively. (d) and (e) The same as (b) and (c) for the error variables $e_{i}$.
Refer to caption: x4.png

## S2 Pseudo-code of MHCACm

The Pseudo-code of Metropolis-Hastings adjusted chaotic amplitude control with momentum (MHCACm) is shown in this section.

Figure: Algorithm 1 Metropolis-Hastings adjusted chaotic amplitude control with momentum

## S3 Boltzmann sampling

### S3.1 Verification of Boltzmann sampling at equilibrium

To test our Markov chain approach, we sample from the Boltzmann distribution of the Ising Hamiltonian $V(\boldsymbol{\sigma})=-\frac{1}{2}\sum_{ij}\omega_{ij}\sigma_{i}\sigma_{j}$. We use the Wishart planted ensemble (WPE) (Hamze et al., ) of size $N=18$ spins to generate the Ising coupling matrix. Using the WPE to generate instances of the Ising Hamiltonian allows us to tune the complexity of the energy landscape and also know the ground state a-priori (Hamze et al., ). Additionally, these instances are known to exhibit properties of NP-Hard problems making the ground state difficult to find by heuristic algorithms (Hamze et al., ). Because we use a small problem size ($N=18$) the exact Boltzmann distribution can be determined by brute force. The sampled distribution obtained by MHCACm is compared to the exact distribution at inverse temperature $\beta=0.07$ by calculating the symmetric KL divergence (see Appendix section S3.2). Note that CACm (without the Metropolis-Hastings step) does not sample from the target Boltzmann distribution (see Fig. S2 a). In contrast, MHCACm exhibits exactly the same distribution of energy as the Boltzmann distribution at thermal equilibrium. This is accomplished while retaining the global exploration properties of CACm, making it a potent hybrid heuristic for identifying lower energy states of non-convex energy functions.

In order to evaluate the sampling characteristics of the constructed Markov chain, we aim at drawing samples from the Boltzmann distribution of the Ising Hamiltonian $V(\boldsymbol{\sigma})=-\frac{1}{2}\sum_{ij}\omega_{ij}\sigma_{i}\sigma_{j}$ for a Wishart planted instance (Hamze et al., ) with a size of $N=18$ spins at an inverse temperature of $\beta$.

Figure: Figure S2: (a) Exact Boltzmann distribution (green dotted line), sampled distribution using CACm (without MH step, see red line), and MHCACm (without MH step, see blue line). Gaussian filtering of the energy distribution is applied for the sake of readability. $\beta=0.07$. (b) KL divergence vs. ratio of inverse temperatures $\kappa=\frac{\tilde{\beta}}{\beta}$. (c) Effective inverse temperature $\tilde{\beta}^{*}$ temperature at which the KL divergence is minimized vs. $\beta$. The states used for the sampled distribution are boolean configurations accepted by the MH criterion.
Refer to caption: x5.png

In Figs. S2 (b) is shown that the KL divergence is minimized for a ratio of inverse temperatures $\kappa=\frac{\tilde{\beta}}{\beta}$. Note that $\tilde{\beta}$ is the effective temperature defined in eq. 2 of the main manuscript. The effective inverse temperature $\tilde{\beta}$ at which the KL divergence is minimized temperature in noted $\tilde{\beta}^{*}$. Figure S2 (c) shows that $\tilde{\beta}^{*}$ is a linear function of $\beta$ within the range of parameters considered. Moreover, the addition of the MH step to MHCACm increases the similarity between $\beta^{*}$ and $\beta$ compared to that of CACm.

### S3.2 Simulations for small N 𝑁 N italic_N

In Fig. S2, the distance to the Boltzmann distribution $P(\boldsymbol{\sigma})=\frac{e^{\beta V(\boldsymbol{\sigma})}}{Z}$ ($Z=\sum_{\boldsymbol{\sigma}}e^{\beta V(\boldsymbol{\sigma})}$) is calculated using the symmetric KL divergence. The symmetric KL divergence $D_{\text{SKL}}(P\parallel Q)$ between the distributions $P$ and $Q$ is defined as follows:

$$ $\displaystyle D_{\text{SKL}}(P\parallel Q)=\frac{1}{2}(D_{\text{KL}}(P \parallel Q)+D_{\text{KL}}(Q\parallel P)),$ (S38) $$

where

$$ $\displaystyle D_{\text{KL}}(P\parallel Q)=\sum_{i}P(i)\log\frac{P(i)}{Q(i)}.$ (S39) $$

In the case of Fig. S2, the convergence of the MCMC is measured by the KL divergence between the exact Boltzmann distribution of energy levels obtained by brute force calculations and the estimated distribution of energy levels obtained by simulation of the Metropolis-Hastings adjusted chaotic discrete space sampling. The former computation is possible because the problems size $N=18$ is sufficiently small for the brute force search method. The instance used for Fig. 2 is a randomly generated instance from the Wishart planted set with hardness parameter $\alpha_{\text{WPE}}=0.8$ (Perera et al., ).

### S3.3 Simulations for large N 𝑁 N italic_N

In order to test the convergence of the MCMC in the case of larger $N$, we can utilize a more heuristic metric. Conventional Gibbs sampling algorithm is run for a sufficiently large number of steps for the chain to converge. In the case of $N=60$ and $\alpha_{\text{WPE}}=0.8$, numerical simulations indicate that the convergence occurs within about 1000 time step of the Gibbs sampling algorithm. To observe the convergence, we limit the characterization of the Boltzmann distribution to its first two moments: the mean energy $<H>$ with $<H>=\sum_{\boldsymbol{\sigma}}H(\boldsymbol{\sigma})P(\boldsymbol{\sigma})$ and standard deviation of the energy $\sqrt{<(H-<H>)^{2}>}$. Figure S3 shows the convergence of the mean and standard deviation of the energy under Gibbs sampling.

Figure: Figure S3: Convergence of the mean and standard deviation of energy distribution of a $N=60$ Wishart planted instance with $\alpha_{\text{WPE}}=0.8$ under Gibbs sampling. Approximate convergence is reached after $1000$ Monte Carlo steps. $\beta=2.5$.
Refer to caption: x6.png

The distance D between the mean and variance of the distribution obtained by the proposed Metropolis-Hastings adjusted chaotic discrete space sampling and the one obtained by iterating the Gibbs sampling algorithm is then calculated as follows:

$$ $\displaystyle D$ $\displaystyle=|<H_{\text{chaos}}>-<H_{\text{Gibbs}}>|$ $\displaystyle+|\sqrt{<(H_{\text{chaos}}-<H_{\text{chaos}}>)^{2}>}-\sqrt{<(H_{ \text{Gibbs}}-<H_{\text{Gibbs}}>)^{2}>}|,$ (S40) $$

where the subscript denote the cases of Gibbs and MH chaotic sampling, respectively.

In Fig. S4 is shown distance $D$ vs. the number of steps $T$ in the algorithm. One step corresponds to one tentative spin flip and one Euler step in the Gibbs sampling and chaotic sampling case, respectively. The proposed sampling method exhibits a faster decrease of $D$ indicating a quicker relaxation to the steady-state distribution. The use of auxiliary variables speed up the convergence compared to the case without.

Figure: Figure S4: Distance D between the mean and variance of the distribution obtained by the proposed Metropolis-Hastings adjusted chaotic discrete space sampling and the one obtained by iterating the Gibbs sampling algorithm (see eq. S38) vs. the number of steps $T$ in the algorithm. One step corresponds to one tentative spin flip and one Euler step in the Gibbs sampling and chaotic sampling case, respectively.
Refer to caption: x7.png

For $N$ larger than $N=60$, measuring the convergence of the distribution to its steady state becomes non-trivial because of strong broken ergodicity (Bernaschi et al., ). Numerical simulation of the convergence to the steady state distribution in this case falls out of the scope of this paper. Note that the convergence to the steady-state distribution can be further accelerated by considering annealing of temperature or the use of the replica exchange Monte Carlo method in both cases of Gibbs sampling and proposed chaotic sampling.

## S4 Autotuning using dynamic anisotropic smoothing

The parameters of the solvers are described in Table 5. We have run the autotuning algorithm called dynamic anisotropic smoothing (DAS) (Reifenstein et al., ) for each problem setting (problem size $N$, bias $b$). The MH sampling rate is set to $f_{\text{MH}}=0.1$.

**Table 5: Parameters of the algorithms tuned using DAS**
| Algorithm | number of parameters | parameter list |
| --- | --- | --- |
| SA | 3 | $\beta_{\text{ini}}$, $\beta_{\text{fin}}$, $T$ |
| AIM | 4 | $\beta$, $\alpha$, $\gamma$, $T$ |
| CAC | 5 | $\beta$, $\alpha$, $\xi$, $a$, $T$ |
| CACm | 6 | $\beta$, $\alpha$, $\xi$, $a$, $\gamma$, $T$ |
| HMCACm | 7 | $\beta$, $\kappa$, $\alpha$, $\xi$, $a$, $\gamma$, $T$ |

An example of the the evolution of parameters during the DAS autotuning is shown in Fig. S5.

Figure: Figure S5: Evolution of parameters during the DAS autotuning for MHCACm. $N=100$, $b=12$, $T=1000$. Here “lamb” denotes $\alpha$ and “kappa” is $\frac{\tilde{\beta}}{\beta}$
Refer to caption: x8.png

## S5 Contruction of degenerate Wishart planted ensemble

The Wishart planted ensemble instances are generated by first choosing a planted ground state $g\in\{-1,1\}^{N}$. However, in our case we consider a set of planted ground states which are represented as columns in an $N$ by $D$ matrix $G\in\{-1,1\}^{N\times D}$. Let $\hat{G}$ be a matrix with the same columns space as $G$ whose columns are orthonormal. Then we generate an $N$ by $M$ matrix $\hat{W}$ of i.i.d. Gaussian random variables. This matrix is then transformed so that the columns are orthogonal to all of the ground states to form the matrix given as follows: $W=(I-\hat{G}\hat{G}^{\top})\hat{W}$. The Ising coupling matrix is constructed as $\hat{J}=WW^{\top}$ with the diagonal elements being set to zero. We can see that every column of $G$ is a ground state since $x^{\top}\hat{J}x\geq 0$ and $G^{\top}\hat{J}G=0$.
To add bias to the energy landscape, we want to somehow create the instance such that configurations close to one ground state tend to have a larger energy than the other ground state even though both ground states have the same energy. To do this, we can consider $W=(I-\hat{G}\hat{G}^{\top})A\hat{W}$ where $A$ is a symmetric matrix and $\hat{W}$ is the random matrix defined above. In this work we consider $A$ to be defined as $A=I+\frac{b}{N}\textbf{1}\textbf{1}^{\top}$ where $b$ is the "bias" away from the ferromagnetic solution. The point of this is that the residual energy (energy away from the ground state) will be amplified for states that are close to the ferromagnetic solution. Then, when choosing the planted ground state we can choose one which is close to the ferromagnetic solution and one that is far creating bias in the search towards one ground state over another. For the numerical results in this work we consider one ground state with a hamming distance of 3 from the ferromagnetic solution and one ground state chosen completely randomly from $\{-1,1\}^{N}$. A Gauge transform is applied to obfuscate the ground states at the end of this process.

## S6 Computing environment

All simulations are done in CPU using a Intel Core i9-11950H, 8 cores, 2.6 GHz and GPU using a NVIDIA RTX A300; except for Table 3 of the main manuscript for which the GPU used is a Nvidia V100 and the CPU wall time is calculated on a 2024 MacBook Pro (Apple M3 Chip). The code is written in python and pytorch.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: NA
soundness: 3
strengths: - The empirical results include multiple baselines and comparisons in terms of different metrics.
- The visualization in Fig.1 clearly shows the main algorithmic idea.
summary: This paper proposes a new algorithm that combines chaotic search and Metropolis-Hastings. The goal seems to solve optimization problems in discrete non-convex energy landscapes. The proposed algorithm is tested on several combinatorial optimization tasks.
weaknesses: - The problem that this paper aims to solve is vague. Is the goal to develop an algorithm that samples better in non-convex energy landscapes, or for optimization in discrete landscapes? Or is the goal to understand non-equilibrium dynamics in non-convex energy landscapes? Similarly, the motivation of the proposed algorithm which combines chaotic search with Metropolis-Hastings is not well-explained.
- The novelty of the proposed algorithm is unclear. Is the algorithm a straightforward combination of chaotic search and MH? If not, what is the challenge, and how does the paper solve the challenge? 
- The empirical improvement is not consistent. For example, Fig.2 shows that CACm is better than proposed method also the variance of the proposed is significantly larger than the baselines.
- The runtime comparison only considers simulated annealing. It will be better to include other baselines as well.
- The paper compared the standard Gibbs with gradient which is developed for combinatorial optimization. It will be more convincing to compare with gradient-based discrete MCMC that is developed for CO, such as [1].

[1] Revisiting sampling for combinatorial optimization, ICML 2023

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: (1)  While the method achieves a 100x speedup over simulated annealing on GPUs, a discussion of any limitations or computational trade-offs encountered in specific scenarios (such as highly multimodal landscapes) would be beneficial, for instance, I wonder if the authors could provide specific examples of problem types or landscapes where their method may face challenges.

(2) I wonder if the authors could provide additional insights into how MHCACm scales with increased problem complexity, for instance, if the authors could demonstrate how the empirical scaling results and the performance of proposed algorithm changes with increasing complexity for a range of benchmark examples.

(3)  The paper implies that the method’s bias towards “easy” ground states is advantageous, but this effect could also limit the algorithm’s ability to reach more challenging or rare ground states. I wonder if the authors could provide quantitative results on the algorithm's performance in finding both "easy" and "hard" ground states across different problem instances. 

(4) Additionally, I wonder if the authors could discuss potential modifications to the algorithm that could help balance the performance of exploration of both easy and hard ground states.
soundness: 3
strengths: (1)The proposed method effectively leverages the Metropolis-Hastings method within a continuous-discrete framework, enhancing sampling efficiency for ground-state discovery in challenging discrete landscapes. This approach demonstrates practical advantages, such as notable improvements in convergence speed and sampling accuracy on GPU architectures.

(2) The method’s focus on non-equilibrium dynamics and its capacity to identify accessible ground states faster than traditional approaches offer a valuable contribution to optimization in rugged energy landscapes.
summary: Overall, this paper proposes a promising hybrid continuous-discrete sampling framework that demonstrates clear benefits in convergence speed and sampling accuracy for rugged energy landscapes. This paper could benefit from additional theoretical insights and comparisons with established methods.
weaknesses: (1) Although the paper demonstrates the practical benefits of the hybrid continuous-discrete approach, the theoretical understanding of the sampling properties of the MHCACm algorithm remains unaddressed. I wonder if the authors could provide a discussion on potential directions for analytical proof of the sampling capabilities of MHCACm, such as convergence rates or mixing times. 

(2) I would like to suggest the authors include more comparison with other prominent sampling algorithms using collective variables, for example, 'Sampling metastable systems using collective variables and Jarzynski–Crooks paths' by G. Stoltz et al. 

In particular, I am curious  to see how the use of collective variables in that work relates to or differs from the proposed method of this paper, and if the authors can combine both approaches.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: * Based on the algorithmic design in this paper, is there any insight we can draw on what an ideal optimizer for deep neural nets should look like?
* Can you elaborate more on the numerical performance of CACm and MHCACm? From the charts, the two performances seem to be close to each other.
soundness: 3
strengths: * The paper writing is nice and structured, with a comprehensive literature review and detailed problem set-up.
* The paper looks technically sound with solid mathematical proof.
* The proposed algorithm is evaluated on multiple tasks and compared with various other benchmark methods, showing competitive performance.
summary: The paper proposes a class of hybrid continuous-discrete algorithms by integrating continuous dynamics with Metropolis-Hastings steps. The paper also constructs a set of Ising problems with a tunable parameter to trade off between easy ground states and hard degenerate ground states, in order to experiment with the bias of different algorithms. The proposed class of algorithms are also fast solvers that achieve a great amount of acceleration on GPU due to a parallelizable structure.
weaknesses: * I am not very familiar with the literature, but seems that the tasks of ground-state sampling are not formally defined in the paper, as well as the idea of non-equilibrium dynamics.
* The connection between ground-state sampling and deep learning optimization/generalization mentioned in the paper is interesting, but the discussion is very limited. 
* For numerical experiments, the definition of TTS is hard to comprehend. Does smaller TTS indicate better algorithmic performance?

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: - The paper only reports the TTS in the experiments on dWPE instances (section 4.4). It would be more transparent if the success probability and runtime data could be provided as well, as it is not clear whether the advantage comes from a higher success probability or a shorter wall-clock runtime due to GPU parallelization. 
- Can you elaborate more on the "dual-primal Lagrangian approach" and its difference from CAC?
- Some minor typos: e.g., line 122 "in order (to) benchmark this algorithm's ability".
soundness: 3
strengths: - State-of-the-art algorithms exploiting relaxation to a continuous state form discrete combinatorial optimization do not sample fairly from the Boltzmann distribution due to the lack of detailed balance. The MHCACm algorithm fills in this conceptual gap by adding a Metropolis-Hasting step. This new design ensures that MHCACm samples fairly from a discrete distribution while iterating over the relaxed (continuous) search space. 
- The numerical results look strong. Table 3 shows that MHCACm has a success probability higher than dSBM, another well-known Ising solver based on GPU. 
- MHCACm is well-suited for large-scale deployment on GPU because its computational bottleneck is matrix-vector multiplication.
summary: This paper presents two variants of the CAC (chaotic amplitude control) algorithm, namely CAC with momentum (CACm) and Metropolis-Hastings CAC with momentum (HMCACm). CACm is a deterministic continuous-time dynamical model for combinatorial optimization, and HMCACm is a Metropolis-Hastings adjusted version of CACm with the Boltzmann distribution as the theoretical equilibrium. MHCACm can be regarded as a unified framework that generalizes many existing methods, including simulated annealing, Hopfield neural networks, analog iterative machines, and CAC(m). Numerical results show that this method illustrates faster relaxation time on NP-hard problems. In particular, MHCACm exhibits excellent performance in sampling from easier ground states, which may be relevant to training over-parametrized neural networks.
weaknesses: - Little theoretical justification for the effectiveness of MHCACm is provided. Specifically, it is not clear why a fair sampling strategy in the CAC framework would lead to a better performance in an optimization (ground state finding) problem. Relaxation to the ground state does not necessarily need to go through a detailed-balance algorithmic path. It would be nice to discuss how the Metropolis-Hastings step interacts with the CAC dynamics to potentially improve optimization performance.
- No empirical results for real-world optimization problems. While the performance of MHCACm has been benchmarked over the dWPE instances and GSET, these test instances are highly artificial and may not reflect the performance of the algorithm in a practical setting (e.g., quadratic assignment problems, portfolio optimization problems, etc.).

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: See Weaknesses.
soundness: 2
strengths: 1. The combination of chaotic dynamics with an MH step to ensure convergence to a target distribution represents an innovative approach to hybrid sampling.
2. The hybrid algorithm shows potential for improving ground-state sampling performance in combinatorial optimization
3. The authors address parallelization for efficient computation on GPUs
4. The discussion on combining chaotic dynamics with probabilistic methods provides a useful context for researchers working at the intersection of machine learning and statistical physics.
summary: The authors propose a novel hybrid continuous-discrete algorithm that combines deterministic continuous dynamics with Metropolis-Hastings (MH) steps for ground-state sampling in non-equilibrium dynamics.
weaknesses: While the authors present an interesting optimization algorithm, the clarity of the writing is a major concern. The main ideas are difficult to follow in the current presentation. I would encourage the authors to reconsider their notation and improve their writing to convey their ideas more effectively to readers.


**Major concerns**
1. The momentum and pre-conditioning typically serve different roles in optimization: momentum accumulates past gradients, and pre-conditioning captures curvature information. I do not think the connection demonstrated in Section 3.2 is trivial. A detailed explanation is needed in the main text to connect them.
2. The paper lacks theoretical guarantees regarding the convergence or performance of the proposed algorithm.
3. Since the algorithm incorporates a momentum variable, it would be more consistent to account for this momentum within the MH step (Equation (7)), rather than applying it solely to $\boldsymbol{\sigma}$.

**Other suggestions**
1. The paper lacks clear definitions for essential variables (e.g., $\boldsymbol{u}$ and $\boldsymbol{e}$ in Equation (1)).
2. The time variable $t$ is used ambiguously. It denotes continuous evolution in Equation (1) but has discrete updates in Equations (4)-(5).
3. Although the authors appear to focus on ground-state sampling, the formulation provided in Section 3.2 is more oriented toward sampling from a Gibbs measure rather than explicitly defining the ground-state sampling.
4. Ensure that all variables and abbreviations are clearly defined. For example, abbreviations such as SA, HNN, and CACm used in Table 1 should be explicitly explained in the main text.

Meta-review: This paper considers designing non-equilibrium dynamics for computing the ground state of rugged energy landscapes. Two variants of the CAC (chaotic amplitude control) algorithm, namely CAC with momentum (CACm) and Metropolis-Hastings CAC with momentum (HMCACm) were proposed. Empirical results show fast convergence of the methods. Reviewers expressed concerns about the precise definition of the problem to solve, comparison with existing approaches, and theoretical justification, but after discussions they suggested most of the concerns were resolved. My impression is overall the strengths overweight the weaknesses, hence the recommendation of acceptance. However, the authors should account for the discussions in a revision.

Decision: Accept (Poster)