Non-Equilibrium Dynamics of Hybrid Continuous-Discrete Ground-State Sampling
Authors: Timothee Leleu
Sam Reifenstein
Venue: ICLR 2025 Poster (2025)
Keywords: Combinatorial optimization
Degenerate ground-state sampling
Metropolis-Hastings algorithm
Chaotic dynamics
Wishart planted ensemble
Primary area: optimization

Abstract: We propose a general framework for a hybrid continuous-discrete algorithm that integrates continuous-time deterministic dynamics with Metropolis-Hastings (MH) steps to combine search dynamics that either preserve or break 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 more quickly, indicating a stronger bias toward these solutions in algorithms using reversible transition probabilities. To validate this, we construct a set of Ising problem instances with a controllable bias in the energy landscape that makes certain degenerate solutions more accessible than others. The constructed hybrid algorithm demonstrates significant improvements in convergence and ground-state sampling accuracy, achieving a 100x speedup on GPU compared to simulated annealing, making it well-suited for large-scale applications.

Published as a conference paper at ICLR 2025
NON- DYNAMICS OF HYBRID
GROUND-STATE SAMPLING
Timothée Guillaume Leleu
Sam Reifenstein NTT Research
NTT Research, Stanford University
t timothee. eeleunnt- research. com t Lleleu@stanford.ees
S samuel rei fensteinnntteeteeeeeeee
ABSTRACT
We propose a general framework for a hybrid continuous-discrete algorithm that integrates continuous-ting deterministic dynamics with Hutropolis-Hastingg (MH) steps to combine search dynamics that either preserve or break detailed balance. Our purpose is to study the non-equilibrium dynamics that leads to the ground state rugged energy landscapes in this general setting. Our results show that MH- driven dynamics reach "easy" ground states more quickly, indicating stronger bias toward these solutions in algorithms using reversible transition probabilities. To validate this, construct set of Ising problem instances with controllable bias in the energy landscape that makes certain degenerate solutions more accessible than others. The constructed hybrid algorithm demonstrates significant improvements in convergence and ground-state sampling accuracy, achieving a 100x speedup on GPU compared to simulated annealing, making well-suited for large-scale applications.
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 al. 22200] Lucas 2014 Mohseni et al., 2022). 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., 2018, /aitt-JJsi et al.] 2018 Feng & Tu [021]: Baldassi et al. 2022, 2023). Here, we study this issue from the general perspective of non-convex optimization, focusing on the Ising Hamiltonian. In particular, we study the task of sampling from degenerate ground states of a non-convex energy landscape. Ground-state sampling involves finding not just any ground state, as often defined in combinatorial optimization, but multiple degenerate ground states. In this context, non-equilibrium dynamics refers to the processes through which systems evolve over time toward steady states, with varying rates of convergence distinct ground states that deviate from equilibrium expectations.
The Metropolis-Hastinnn (MH) algorithm is one of the most widely used techniques for sampling from reversible Markov Chain Monte Carlo (MCMC) (Zhang al.) 2020 Chib & Greenberg 1995 Chen et al. 2014). This approach proposes new states and accepts them with 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 these energy landscapes often arise in systems exhibiting broken ergodicity (Bernaschi et al. 2020). such as spin glasses. In these cases. other MCMC algorithms that maintain detailed balance, like simulated annealing (Kirkpatrick et al. 1983) and more advanced methods (Hukushima & Nemoto 1996) Surraaovvic al.) 2022 Hukushima & Iba 2003 Houdayer] 22001] Landau et | 2004). 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. 2019 2021 Goto et al. 2021; Buesing et al. 2011). particularly those based chaotic search (Leleu al. 2019 2021). Although these approaches do not sample in general from the exact Boltzmann distribution at equilibrium, can be shown in special cases that non-reversible Markov chains mix faster than reversible ones (Kapfer & Krauth

---

Published as a conference paper at ICLR 2025
2017). However, their ability to reach degenerate ground states of a rugged energy landscape using non-equilibrium dynamics is still not well understood.
In this work, we introduce generalized algorithm that combines MCMC techniques with chaotic search dynamics via the chaotic amplitude control (CAC) algorithm. Tuning the hyperparameters of this hybrid approach allow for smooth transition between the two strategies. Additionally, create new planted instances with degenerate ground states, inspired by the Wishart planted models. 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 's performance can be optimized based on the bias parameter, making it competitive with state-of oftte-of-the-ar methods. The code for our algorithm can be found on this online repository: GitHub Repository
RELATED WORK
2.1 DISCRETE STATE MCMC OPTIMIZERS
MCMC methods for optimizing energy functions over discrete spaces-such as simulated annealing
(Kirkpatrick et al.] 1983) and parallel tempering (Hukushima & Nemoto 1996) are standard approaches in the field of optimization. Building upon these, recent schemes such as population annealing (Hukushima & Iba 2003), variational parallel tempering (Surjanovic et al. 2022) and specialized techniques, including cluster updates (Houdayer) 200T), Wang-Landau sampling (Landau ctt al. 2004), and Replica Exchange Wang-Landau Sampling (Vogel et al. (0013), 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. This is a behavior observed in many systems with frustrated interactions. (Bernaschi et al., 2020).
2.2 HYBRID CONTINUOUS-DISCRETT 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 1998), Hamiltonian Monte Carlo (HMC) (Neal et al. 2011; Xifara et al. 2014), as well as derivative methods such as the No-U-Turn Sampler (NUTS) (Hoffman et al. 2014). Generalized HMC (Chen
et al., 2014). HMC with the Gaussian integral trick (Zhang et al. 2012), amortized Metropolis adjustments (Zhang et al., 2020). and sampling using collective variables (Schönle et al. 2024). These methods primarily target continuous spaces, and thus are not directly applicable 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. 2020). Gibbs with Gradient (GWG) (Grathwohl et al. (02T), and the Hamming Ball Sampler (Titsias & Yau 2017).
Nonetheless, these approaches not suitable for ground-state sampling.
2.3 HYBRID
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 al., 2002)), all designed to solve optimization problems more efficiently than current nigh-perffrmmnce digital computers.
In particular, dynamical systems such as simulated bifurcation machines (Goto et al. 2021 analog iterative machines (AIM) (Kalinin et al.) (023), and chaotic amplitude control (CAC) (Leleu et al.) 2019 |021 Reifenstein et al.) (023)) offer an alternative approach. Rather than relying on gradient descent of an energy function in systems with symmetric interactions, such in the seminal work by Hopfield (Hopfield & Tank 1985). it has been shown that Hamiltonian dynamics (Goto et al. 2021 Neal et al. 2011) or chaotic dynamics driven by asymmetric interactions (Leleu al.) 2019 2021) can outperform standard heuristics like simulated annealing (Kirkpatrick et al. 1983), parallel
tempering (Hukushima & Iba. 2003), and breakout local search (Benlic & Hao 2013) on specific
2

---

Published as a conference paper at ICLR 2025
benchmark datasets (Leleu al., 2021). This has challenged 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 transi- tions 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 1970) (see Fig. [a). The algorithm built on three main concepts: (1) relaxation from discrete to continuous space (or continuous embedding), (2) addition of auxiliary variables to es- cape from local minima (chaotic amplitude control (Leleu et al. 2019)) 202I) and momentum (Kalinin al. 2023 (3) probabilistic jumps in the discrete space based the Metropolis- -asstings criterion.
In the next section, we construct a novel set of Ising problem instances inspired from the Wishart planted ensemble (Hamze al.) 2020 in order benchmark this algorithm's 's ability 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. [b).
b
Chaotic Netropolis-Hastingg jump
More likely sampled
Lesa likey sempled
Energy landscagpe
Ground-state anergy
Continuous
Discreti configuration space
Figure 1: a) Proposed hybrid approach generates deterministic chaotic trajectories that rapidly explores the energy landscape and probabilistic jumps accepted via the Metropolis Netropolis-Hasttnns criterion. (b) Constructed planted instances with degenerate ground states with unequal sampling probabilities by heuristics. Some ground states are easier, and others harder, for local search heuristics to discover.
3.2 DETERMINISTIC CONTINUOUS TIME MODEL
We define a continuous time dynamical system inspired by the CAC algorithm (Leleu al. 2019 2021) and generalize its formulation to draw a correspondence between its parameters and the inverse temperature 3 of the Boltzmann distribution P(ar) x eV()), The goal is to design a dynamical system capable of sampling from the ground states of V. specifically from the sero-temperature or spins.
distribution P(ar) defined on the discrete space o E {-1, 13 where N is the number of variables
First, discrete variables are relaxed to the continuous space E V and we define the potential (or energy function) V(x) on the real space. The following dynamical system is used direct the search towards lower energy states:
duu
du di2
-cu - eo V.V,
di
where the symbol denote the Hadamard product; V, V. the gradient of V with respect to the vector æ; a, positive parameter, and e, a vector of positive auxiliary variables. The variables are often referred to "soft spins." The variables represent the internal states of these soft spins, while the auxiliary variable accounts for variations in their amplitude (Leleu et al., 2019). The term proportional to y represents the momentum which has been utilized for non-convex
3

---

Published as a conference paper at ICLR 2025
optimization (Kalinin et al. 2023) and sampling (Neal al. 201 Xifara et al. 2014). Moreover, we set:
24 ((Pen))), Vi E ,
and ((x)) - I+e- is sigmoidal function normalized to the domain x E [-1, 1]. Serr is an effective temperature (see Appendix section S3.1 for a discussion about its influence).
2
The auxiliary variables e; undergo dynamic changes over time, with the goal of rectifying amplitude heterogeneity, given as follows (Leleu et al. 2019):
de dt
E(x: o a a) oe.
with and being parameters: represents target amplitude for the variables x and & represents the speed of the error correction dynamics. Because this algorithm adds momentum dynamics to the original chaotic amplitude control equations (Leleu al., 2019 2021). we call chaotic amplitude control with momentum (CACm).
In the Appendix section SI, the role the auxiliary variables e is analyzed and we show that their modulation can be interpreted as dynamic pre-conditioning of the Hessian of V at the proximity of the discrete state. Because the auxiliary variables 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 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. Note that the momentum and pre-conditioning play distinct roles in the dynamics. The momentum aggregates past gradients, while pre-conditioning modulates the effective curvature of V.
3.3 PROBABILISTIC DISCRETE TIME MODEL AND MCMC
In order to utilize this dynamical system algorithmically, discretize eqs. 2 using straightforward Euler approximation. Unlike Hamiltonian Monte Carlo, which requires careful selection of the integration method, accurate integration scheme is not essential in this context since the goal is to sample from a discrete distribution rather than precisely from a continuous one. 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 defined follows (after suitable change of variables):
u(m + 1) = u (m) - am((m) eV eV (m(m)) + (u((m) u(m - 1)), e(m 1) e(m) {(x((m) x (m)
(4)
e(m).
(5)
To maintain the effect of the deterministic dynamics when used for optimization, we formulate an MCMC as follows. The eqs. and [] iterated for n time steps, subsequent to which we sample discrete state o f-11, 1) as detailed below:
P(o, = T +
0, = otherwise.
(6)
2
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 a can leap remote states. Note moreover that when we set & 0, e0) = 0, 1, = 1, and update spin at time, we revert to the equation of the Boltzmann machine.
3.4 METROPOLIS ADJUSTED CHAOTIC SAMPLING
When & / 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 in principle) to the Boltzmann distribution, we introduce samples okk)) where & belongs to {1, 21., K}. The inverse of the number of Sattroools--Hssting steps is denoted SMu k and the total number of steps is noted T' with T Kn.
dietropoiis-Hastine
acceptance criterion for the discrete state
4

---

Published as a conference paper at ICLR 2025
Typically, the acceptance criterion of new configuration т, when the current state is is given as follows:
P(T)Q(O/T) P(T))(
A(T/a) min(1,
Here, Q(T/or) represents the probability of transitioning from to The Monte Carlo step for our
method is comprised of two substeps as follows: (1) A deterministic trajectory following eqs. f Mlad
for n steps starting from u(kn) = 0. ekkn) = Õ, and <<kn) = o (kn)) (2) A sample T generated
with the probability p with p y+1 with defined as ys = ø( Senu, Sen,, ((=+()))))))))))) ), i.e., P(r; = = ps. ae ((k+(k)))
-1 otherwise.
Note that y depends exclusively on o, since it is generated deterministically from it, Because of the synchronous updating of all spins, we can write the following definition for Q(T/or):
1 Ti log (1 - pi) 2 To log (p:) +
log Q(T/O)= = Σ
(8)
2
To calculate the transition probability in the reverse Monte Carlo step from T to σ, we need to apply
first the deterministic sub-step (1) from to obtain candidate end point af from which we can calculate the probability to generate the random sample o with probability p #+1 The acceptance
2
criterion can then be rewritten follows:
N
log A(T/o) = min(0, min(0,,(((o))))) 3(V(o) - V(T)
Σ
+ Ti
log (p:) +
1 Ti log(1 (1 - p.)
2
2
1+0 o log (1)
log (1 pi)
2
This acceptance criterion can be computed efficiently since it only requires the logarithm of terms already computed within the deterministic equations. The discrete state T is accepted with probability A(T/o). The Markov chain using the Metropolis-Hastings step respects detailed balance condition and converges to our target Boltzmann distribution P(x) (z)
In summary, the pscudo code of our approach is given in section S2 the Appendix. Note that y is solely determined by The deterministic path to y originates from with 1. Consequently, the transition probability from a to T depends only on a (mapped to y 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. | We call the algorithm defined by eqs. 49 Metropolis Hastings Chaotic Amplitude Control with momentum (MHCACm).
3.5 SUMMARY OF ALGORITHMS
key property of new the algorithm proposed in this work is that can be interpreted generalization of many existing methods which exist as limiting cases of MHCACm. Some examples are: Simulated annealing (Kirkpatrick et al., 1983). Hopfield neural networks (Hopfield & Tank 1985). analog iterative machines (Kalinin et al. 2002), chaotic amplitude control (Leleu et al. 2019), ,, and chaotic amplitude control with momentum. These limiting cases are summarized in Table 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 performs at least well as any of its constituent parts. Importantly, automatic parameter tuning methods have recently been developed which allow the many parameters to be quickly and accurately tuned (Reifenstein et al. 2024).
defined similarly to y using the end state of the deterministic path starting from given as 2; Serts ((b+2)nn) exe, ((k+2)))
5

---

Published as a conference paper at ICLR 2025
Table 1: Limit cases of the proposed algorithm MHCACm. SA: simulated annealing. HNN: Hopfield
neural network. AIM: analog iterative machine. CAC: chaotic amplitude control. CACm: chaotic amplitude control with momentum. MHCACm: Metropolis Hastings Chaotic Amplitude Control with momentum.
model
fwn
E
model
Y
SA
= 7
=0
CAC
= 0
> 0
HNN
= 1
=0
CACm
1
> 0
0
AIM
=1
>
MИСАCm
<
Sam <1
> 0
> 0
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. 2020). WPE instances are particularly useful for benchmarking due to their tunable hardness parameter OswPE- For certain values of owpe: the recovery of the planted solution becomes easier or harder. We introduce the degenerate Wishart planted ensemble (dWPE) to allow for planting of multiple degenerate ground states. Additionally, we add 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. (b). Thus, we refer these different ground states as "easy" and "hard" solutions. Details of the construction method are provided 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 combina- torial 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 (the smaller, the better the algorithm's 's performance). 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 ground state before fully equilibrating to the Boltzmann distribution and, consequently, each degenerate solution can have different probability of being sampled at finite sampling time. We define TTS, TTSaay, and TTSpad as the time to reach any, the easier, and harder find degenerate solution (of same Ising energy) with 99% probability as follows:
TTS = log(1 1 0.99)
log (1 Σ, TTS easy = T log(1 (1 0.99)
log (1 max.Pa) TTS = T log(1 0.99) (10) log (1 min.Pa)
where Pa is the probability to sample the ground state d. The algorithms used for comparison are all tuned to nearly optimal parameters by minimizing the TTS. We use a state-of-the-art autotuning method called dynamic anisotropic smoothing (Reifenstein et al.) (024) (DAS) which is particularly well fitted to optimizing the parameters of heuristics for combinatorial optimization (Reifensteim et al. 2024). 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 single step is dominated by the matrix-vector multiplication of size N N
4.3 GLOBAL OPTIMIZATION TO ANY OPTIMAL SOLUTION
For the task of finding any ground state, Fig. (a) shows that MHCACm performs almost well as CACm and better than AIM, CAC, and SA In the case of sampling any ground state of biased
For SA, a step represents full sweep of N spin updates which also has complexity O(N2). The for PT, but divided by the number replica.
Results for SA and PT are collected using dwave-neal (dwavesystems 2024 and PySA Mandra & Munoz
Bauza 2024

---

Published as a conference paper at ICLR 2025
(b = 12) dWPE problem instances (see section 4.1 than other algorithm to find an optimal solution (see Fig. (b)).
MHCACm requires significantly fewer steps
(unbiased Wishart)
b (biused Wishart b 12)
107
10
10
10F
TTS pryggggnngnn
10
10
10"
10*
TTS
107
CAC
102
10a
60
100
130
120
130
100
110
1.20
150
N
Figure 2: The time to solution of (unbiased) Wishart planted instances and biased degenerate planted instance with b = 12 are shown (a) and (b), respectively. fau 0.1. All other parameters tuned with DAS (Reifenstein et al. (024) (see Appendix S4). Colored regions represent the 95% confidence interval calculated by generating 15 instances per data point.
4.4 SAMPLING IN THE CASE OF BIASED DEGENERATE OPTIMAL SOLUTIONS
Next, examine why MHCACm requires fewer steps to find any solution when there is larger bias in the dWPE instances. Fig. (a) and (b) show that the TTS for easy ground states decreases function of the bias b of the generated Wishart instances. Conversely, the time sampled the hard solution TTStarrd increases. This is unlike the case of CACm without the MH step, for which the TTS and TTS 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 Iindicates 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 The energy landscape contains higher density of low-energy states near the easy solution. 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 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.
We show the results of the Boltzmann sampling algorithm Gibbs with gradient (Grathwohl et al. 2021) (GWG), although its purpose is not combinatorial optimization. In appendix section S6. comparison another recently proposed sampling algorithm (Sun et al. (003) is shown.
We also investigate the effect of changing the complexity parameter OwPF of biased Wishart planted instances (see Fig. with bias of = 12. We observe that the relative reduction of the TTS due to the introduction of the MH step in biased instances is more pronounced for instances of higher complexity (i.c., smaller parameter awPE as shown by the comparison of Fig. (a) and (c)).
4.5 RESULTS ON GSET
The GSET Ye 2024 is set of Max- Cut problem instance which are commonly used to benchmark
Ising and Max- Cut solvers (Goto et al. 2021 Leleu et al. 2019 Benlic & Hao 2013). To test the ability of MHCACm to find optimal solutions on larger optimization problems, we briefly present some results these instances in comparison with state of the art algorithm known as dSBM (Goto
7

---

Published as a conference paper at ICLR 2025
TTS, "(ten) grunn-uoon
109
0.5
175ea (63ty gound
10
10
1000
2000
3000
4000
101
101
107
TTSaa
10"
10
CACH
-K- CAC
10
CACa:
-K-
CAC
MERCAOE
-4: P7
MEROAOK
ADN
AIN
30"
10
1000
2000
4000
bsech
Figure 3: Time to reach the easy solution TaSy (a) and hard solution TTS (b) vs. the bias of degenerate Wishart planted instances. (c) max runtime T for degenerate Wishart planted instances. The same optimal parameters as in Fig. are used. (d) The same for TTSa (a,b,c,d) N 100. (c,d) b 10.
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, cay and T (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, TTSd not defined for it and hard TTS. Bold numbers denote the best algorithm.
10
TTS
10 3 TTSaus
10
TTSaa
model/N
100
140
100
140
100
140
SA
27.56
1917.67
27.56
1917.67
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
et al.
2021 1
Table 3 shows that MHCACm reaches best known solutions faster than dSBM (Goto on the first few instances of the GSET problem set Ye 2024).
et al. 2021
Table 3: TTS (in steps or MVM) for solving GSET
Ye
2024) instances of MHCACm and
dSBM Goto et al. 2021).
instance
MHCACm
dSBM
N (id)
Pu
TTS
Po
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.1 16500
3000
75858
0.019060
1500
358947
4.6 PARALLELIZATION ON GPU
The algorithm 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 the matrix vector multiplication (MVM), a basic linear algebra operation that can be done very efficiently modern GPUs using frameworks such PyTorch. In table we compare the wall-clock TTS of our PyTorch implementation of MHCACm finding ground states a degenerate WPE instance of size N 100. We also show the wall-clock time of standard implementation of SA (Tiosue
8

---

Published as a conference paper at ICLR 2025
CACm MHCACD
CACm MBCACm
CACm MHCACu
10-TTSS
TIS
10 TTS
10
0.73
53D
0.00
6.50
0.90
0.00
0.70
0.50
0.50
Wishar complexity mane
Wishart complexity mare
Widarrt Compplccool
Figure 4:
TTS, ETS and Thd (as defined in the text) vs. complexity of Wishart planted
instances OWNE (see definition in Appendix section S5). N 100, b = 12.
2024 rumning CPU. MHCACm on GPU exhibits 100x speed up in real time against simulated annealing on CPU for sampling easy optimal solutions. Furthermore, MHCACm is compatible with an implementation special-purpose analog hardware designed for linear algehra acceleration, such as a memresistor crossbar (Hu et al. 2014) or on-chip photonic mesh arrays (Bogaerts et al.) 2020).
Table 4: Wall clock times of algorithms different computing platforms. TTS is time to find any ground state of a degenerate WPE instance with OWPE = 0.8, N 100 and bias b = 12. The GPU used is Nvidia V100 and the CPU wall time calculated 2024 MacBook Pro (Apple M3 Chip).
TTS (s) (N 100)
unbiased WPE
biased WPE (b 12)
SA (CPU)
0.0756
0.0567
AIM (CPU)
0.0287
CACm (CPU)
0.013
0.016
MHCACm (CPU)
0.0163
4.08 x 10
MHCACm (GPU)
0.00207
5.18 x 10
COMPARISON TO OTHER APPROACHES
Methods such HMC (Neal al. 2011 Xifara et al. 2014 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, making them ill-suited for discrete state sampling. CAC suggests different mechanism for the auxiliary variables. adjusts the amplitude of the continuous variables, ensuring proximity to the discrete (binary) state (Leleu et al. 2019 2021). Conventionally, this can be interpreted as a dual-pprmal Lagrangian problem in which we want to minimize the error function (in the continuous embedding) while also satisfying constraint (being close to the discrete state). Our method differs from the primal dual /Vadlamaní et al., 2020) approach. With CAC, the error function is achieved by multiplying the gradient with auxiliary variables, akin to dynamically pre-conditioning the Hessian the cost function.
Unlike HMC however, CAC does not conserve (even approximately) energy. In HMC (Neal et al. 2011 Xifara et al., 2014), the MH step compensates for the discretization errors brought about by the leapfrog integration. However, the reason for adding an MH step CAC is different. State-of-the-an algorithms that exploit relaxation to an continuous state for discrete combinatorial optimization such memcomputing (Sheldon al. 2019) chaotic amplitude control (Leleu al., 2019), and simulated bifurcation machine (Goto et al. 2021) do not accurately sample 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 hybrid MCMC.
6 CONCLUSION
The MHCACm algorithm effectively combines chaotic dynamics with the Metropolis-Hastings method to enhance ground-state sampling efficiency discrete spaces. The introduction of the MH step increases the probability of finding more accessible ground states during non-equilibrium dynamics when a bias exists in the energy landscape. However, the MH step does not seem beneficial when such bias among degenerate ground states is absent. This approach significantly improves the

---

Published as a conference paper at ICLR 2025
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 Appendix Fig. S2) but can exhibit faster relaxation times on NP-hard Ising problem instances with rugged energy landscape (see Fig. and Tab. I when compared 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 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 al. 2018: Baity-Jesi al. 2018 Feng & Tu 2021 Baldassi et al. 2022 2023).
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Marco Baity-Jesi, Levent Sagun, Mario Geiger, Stefano Spigler, Gérard Ben Arous, Chiara Cam- marota, Yann LeCun, Matthieu Wyart, and Giulio Biroli. Comparing dynamics: Deep neural networks versus glassy systems. In International Conference on Machine Learning, pp. 314-323,
PMLR, 2018.
Carlo Baldassi, Clarissa Lauditi, Enrico M Malatesta, Rosalba Pacelli, Gabriele Perugini, and Riccardo Zecchina. Learning through atypical phase transitions in overparameterized neural networks. Physical Review E,
2022.
Carlo Baldassi, Enrico M Malatesta, Gabriele Perugini, and Riccardo Zecchina. Typical and atypical solutions in nonconvex neural networks with discrete and continuous weights. Physical Review E, 108(2):024310, 2023.
Una Benlic and Jin-Kao Hao. Breakout local search for the max-cutproblem. Engineering Appli- cations of Artificial Intelligence, 26(3): 162-1173, 2013. ISSN 0952-1976. doi: https://oi.orr 001. URL httts:: //wp://wwwwwww www.:science... rect. com/ science/ arricleess11/09922177772000
Massimo Bernaschi, Alain Billoire, Andrea Maiorano, Giorgio Parisi, and Federico Ricci-Tersenghi. Strong ergodicity breaking in aging of mean-field spin glasses. Proceedings of the National Academy of Sciences, 117(30):15522-177)
2020.
Wim Bogaerts, Daniel Pérez, José Capmany, David AB Miller, Joyce Poon, Dirk Englund, Francesco Morichetti, and Andrea Melloni. Programmable photonic circuits. Nature, 586(7828)2007-216, 2020.
Lars Buesing, Johannes Bill, Bernhard Nessler, and Wolfgang Maass. Neural dynamics sampling: model for stochastic computation in recurrent networks of spiking neurons. PLoS computational biology, ((11)=e1022211 2011.
Tianqi Chen, Emily Fox, and Carlos Guestrin. Stochastic gradient hamiltonian monte carlo. In International conference on machine learning, pp. 1683-1691, PMLR, 2014.
Siddhartha Chib and Edward Greenberg. Understanding the merropolis-hastings algorithm. The american statistician, 4944)3327-335, 1995.
dwavesystems. dwave-neal. http://gghabbooommmowsssss////ee 2024. Accessed: 2024-08-
15.
Yu Feng and Yuhai Tu. The inverse variance-flatness relation in stochastic gradient descent is critical for finding flat minima. Proceedings of the National Academy of Sciences, 118()))20156671118 2021.
Hayato Goto, Kotaro Endo, Masaru Suzuki, Yoshisato Sakai, Taro Kanao, Yohei Hamakawa, Ryo Hidaka, Masaya Yamasaki, and Kosuke Tatsumura. High-perfone gegh-performnnce combinatorial optimization based classical mechanics. Science Advances, 7()):abbe9533, 2021.
Will Grathwohl, Kevin Swersky, Milad Hashemi, David Duvenaud, and Chris Maddison. Oops took a gradient: Scalable sampling for discrete distributions. International Conference on Machine Learning, pp. 3831-3841 PMLR, 2021.
Firas Hamze, Jack Raymond, Christopher A Pattison, Katja Biswas, and Helmut G Katzgraber. Wishart planted ensemble: A tunably rugged pairwise ising model with a first-order phase transition. Physical Review E. 101(5):05210,, 2020.
W Keith Hastings. Monte carlo sampling methods using markov chains and their applications. 1970.
Matthew D Hoffmman, Andrew Gelman, al. The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo. J. Mach. Learn. Res., 15(1):1593- 2014.
John J Hopfield and David W Tank. "neural" computation of decisions in optimization problems. Biological cybernetics, 52(3):141-152 1985.
11

---

Published as a conference paper at ICLR 2025
Jérôme Houdayer. A cluster monte carlo algorithm for 2-dimensional spin glasses. The European Physical Journal B-Condensed Matter and Complex Systems, 22:479-44 4844,4842220 2001.
Miao Hu, Hai Li, Yiran Chen, Qing Wu, Garrett S Rose, and Richard W Linderman. Memristor crossbar-based neuromorphic computing system: case study. IEEE transactions on neural networks and learning systems, 1811))8864-1875 2014.
Koji Hukushima and Yukito Iba. Population annealing and its application to a spin glass. In A/P Conference Proceedings, volume 690, pp. 200-206. American Institute of Physics, 2003.
Koji Hukushima and Koji Nemoto. Exchange monte carlo method and application to spin glass simulations. Journal of the Physical Society Japan, 65(6): 1546):16041108, 1996.
Kirill P Kalinin, George Mburgias-Alexandrr Hitesh Ballani, Natalia G Berloff, James H Clegg, Daniel Cletheroe, Christos Gkantsidis, Istvan Haller, Vassily Lyutsarev, Francesca Parmigiani, et al. Analog iterative machine (aim): using light solve quadratic optimization problems with mixed variables. arXiv preprint arXiv:2304. 12594, 2023.
Sebastian C Kapfer and Werner Krauth, Irreversible local markov chains with rapid convergence towards equilibrium. Physical review letters, 1 119((4))240033, 2017.
Scott Kirkpatrick, C Daniel Gelatt Jr, and Mario P Vecchi. Optimization by simulated annealing. science, 220(4998)6611668) 1983.
DP Landau, Shan-Ho Tsai, and M Exler. A new approach to monte carlo simulations in statistical physics: Wang-landau sampling. American Journal of Physics, 72(10):1294 199-11322222 2004.
Timothée Leleu, Yoshihisa Yamamoto, Peter McMahon, and Kazuyuki Aihara. Destabilization of local minima in analog spin systems by correction of amplitude heterogeneity, Physical review letters, 122(4))040607, 2019.
Timothée Leleu, Farad Khoyratee, Timothée Levi, Ryan Hamerly, Takashi Kohno, and Kazuyuki Aihara. Scaling advantage chaotic amplitude control for high-performance combinatorial optimization. Communications Physics, 4(1):266, 2021.
Andrew Lucas. Ising formulations of many np problems. Frontiers in physics, 2:74887, 2014.
Salvatore Mandrå and Munoz-Bauza. Pysa: Fast simulated annealing in native python.
https//gitttb.coomnns////// 2024. Accessed: 2024-08-15.
Naeimeh Mohseni, Peter L McMahon, and Tim Bymees. Ising machines as hardware solvers of combinatorial optimization problems. Nature Reviews Physics, 4(6):363-379 2022.
Radford M Neal et al. Mcme using hamiltonian dynamics. Handbook of markov chain monte carlo, 2(1 1):2, 2011.
Akihiko Nishimura, David Dunson, and Jianfeng Lu. Discontinuous hamiltonian monte carlo for discrete parameters and discontinuous likelihouds. Biometrika, 107(2):365-380 2020.
Dilina Perera, Inimfon Akpabio, Firas Hamze, Salvatore Mandra, Nathan Rose, Maliheh Aramon, and Helmut G Katzgraber. Chook-a comprehensive suite for generating binary optimization problems with planted solutions. arXiv preprint arXiv:2005. 14344, 2020.
Sam Reifenstein, Timothee Leleu, Timothy McKenna, Mare Jankowski, Myoung Gyun Suh, Edwin Ng. Farad Khoyratee, Zoltan Toroczkai, and Yoshihisa Yamamoto. Coherent solvers: tutorial. Advances in Optics and Photonics, ((2):335-441 2023.
Sam Reifenstein, Timothee Leleu, and Yoshihisa Yamamoto. Dynamic anisotropic smoothing for noisy derivative-free optimization. arXiv preprint arXiv:2405.01731 2024.
Gareth O Roberts and Jeffrey S Rosenthal. Optimal scaling of discrete approximations to langevin diffusions. Journal of the Royal Statistical Society: Series (Statistical Methodology). 60(1): 255-268, 1998.
12

---

Published as a conference paper at ICLR 2025
Christoph Schönle, Marylou Gabrié, Tony Lelièvre, and Gabriel Stoltz. Sampling metastable systems using collective variables and jarzynski-crooks paths. arXiv preprint 2rXiv2240S.88160 2024.
Forrest Sheldon, Fabio Traversa, and Massimiliano Di Ventra. Taming a nonconvex landscape with dynamical long-range order: Memcomputing ising benchmarks. Physical Review E, 100(5): 053311, 2019.
Daniel Soudry, Elad Hoffer, Mor Shpigel Nacson, Suriya Gunasekar, and Nathan Srebro. The implicit bias of gradient descent on separable data. Journal of Machine Learning Research, 19(70):1-57, 2018.
Haoran Sun, Katayoon Goshvadi, Azade Nova, Dale Schuurmans, and Hanjun Dai. Revisiting sampling for combinatorial optimization. International Conference on Machine Learning, pp- PMLR, 2023.
Nikola Surjanovic, Saifuddin Syed, Alexandre Bouchard-Côté, and Trevor Campbell. Parallel tempering with a variational reference. arXiv preprint urXiv:2206.0000 2022.
Joshua Tiosue. Qubovert: A python library for representing and solving qubo problems. https//ggthhbccomoooueeeUUUUCC 2024. Accessed: 2024-05-22.
Michalis K Titsias and Christopher Yau. The hamming ball sampler. Journal of the American Statistical Association, 1 112(520):1598-16)) 1611, 2017.
Sri Krishna Vadlamani, Tianyao Patrick Xiao, and Eli Yablonovitch. Physics successfully implements lagrange multiplier optimization. Proceedings of the National Academy of Sciences, 117(43):
2020.
Thomas Vogel, Ying Wai Li, Thomas Wüst, and David P Landau. Generic, hierarchical framework for massively parallel wang-landau sampling. Physical review letters, 110(2))211003, 2013.
Tatiana Xifara, Chris Sherlock, Samuel Livingstone, Simon Byrne, and Mark Girolami. Langevin diffusions and the metropolii-adjusted langevin algorithm. Statisties & Probability Letters, 91: 19, 2014.
Yinyu Ye. Gset: A collection of graphs, 2024. URL URLLttpps//mbb..tttodddddd////// Tyye/Gset/ Accessed: 2024-05-22.
Ruqi Zhang, A Feder Cooper, and Christopher De Sa. Amagold: Amortized metropolis adjustment for efficient stochastic gradient meme. International Conference on Artificial Intelligence and Statistics, pp. 2142-2152 PMLR, 2020.
Yichuan Zhang, Zoubin Ghahramani, Amos J Storkey, and Charles Sutton. Continuous relaxations for discrete hamiltonian monte carlo. Advances in Neural Information Processing Systems, 25, 2012.
13

---

Published as a conference paper at ICLR 2025
APPENDIX
EOU-EQUILI DYNAMICS OF HYBRID CONTINUOUS-DISCRETT BROUND-STATE SAMPLING
SI 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 7 =
S1.1 DYNAMIC TRE-CONDITIONING OF THE HESSIAN OF
The equation (1) of the main manuscript is 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 e and explain that their modulation can be interpreted as dynamic pre-conditioning of the Hessian of V at the proximity of the discrete state a. First, we show Appendix section S1.2 that the variables u in the system of eqs. (1) and (2) can be eliminated by change of variables. Then, the Jacobian matrix Jaw in the space æ defined as J.z = { Au; dzy hif with ē; de
di can be expressed as follows at proximity of the fixed points of the
system (see Appendix section SL.3
J = DI
Die g(vor) (ao)
(S1)
where D[x=] is a diagonal matrix with elements given by the components of the vector æ and H(Vaor is the Hessian of V defined as the Jacobian of V and calculated at the point vaor The functions I and and defined in Appendix section S1.22 Equation [ST shows that the auxiliary variables e modulates the Hessian H on the slow time scale such that (t, = D(e(r)M(x;)) with h(x)
for æ E] - 1, 1[ and (xx))
(x) (see lemma in Appendix section $1.3
Ban
Given that the fixed points of eqs. (1) and (2) are given as e; -avaa 2Z (vaos) we can interpret
the role of normalization of each row of the Hessian by the inverse of the gradient in the corresponding direction (see theorem in Appendix section [ST.3). Because the auxiliary variables remain positive, they do not change the inertia (number of positive, negative, and zero eigenvalues) of the Hessian (see lemma in Appendix section (1.44) but can change its directions and eigenvalues. Moreover, 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 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 tftte-of-th-arr heuristics to solve Ising problems (Leleu et al. 2019 2021
$1.2 CHANGE OF VARIABLE
Consider x; = ((2Buu)) = u, = (xx) with (x;;) = (-1(2x)
29
In the case (xx) (x)())) Ova(z) Moreover, (z) is positive on | 1, 1[ and
2 I+e-s 1, we have o (xx) = -In( z
1), (/((x) = Ma+1|a-11 with
for |x| +
Because dus (x,)) des with vi/3 = 1 B0 (a,) eliminating u, (Vi E f1,..., ,N) as follows:
the equations of motion can be rewritten by
di
az.
dr;
(g(z1) C,
V.V.
(S2)
dt V/s(x;) (x(),)
5(0? - a)e.. de,
(S3)
dt
14

---

Published as a conference paper at ICLR 2025
$1.3 JACOBIAN MATRIX THE FIXED POINTS
The fixed points deg as follows:
and
at the dynamical system described by eqs. S2 and are given
at
24 = 0, o.va.
(S4)
(
4 = av Ow: (ava)
(S5)
Moreover, the Jacobian matrix the fixed points is defined as:
Jya
Jue
J
(S6)
Jex
Jee
with
Jan = D[-aa(Vaoa))) - Die [ Dle ((Vao))
(S7)
Jue = D11 (Vaoi)
JV
}:),
(S8)
2 a
Jex 2EvaDlo
(S9)
Jee = 0,
(S10)
with H(x) = V is the Hessian of V in x without its diagonal elements. The notation D(æ)
denote the diagonal matrix with diagonal elements given by the components of the vector The functions f, 9, and of eq. S7 are defined as follows:
Va
(ar)
f(x;)
(S11)
Ox
8V
(x:) Bi:
g(x:)
(S12)
Ox:
1
1(x;)) (x)():))
for H, E] 1,
(S13)
In the case of 6(z)
have:
Ite
1 Di
f(x;) = 2 + x,In(
(S14)
1 + Mi
0V
02V
g(x,) = 3x, dx:
B(xx + 1)(xx - 1)
(S15)
dr?"
Mxx)) = - -((r + 1) - 11)
(S16)
Note that f is an even function with f(x) f(-x).
In the case { + O, the dynamics of e, variables becomes much slower than that of z;. We note r the slower time scale of e; with tE. The lemma is obtained using eq. [S7]
Lemma When ş + 0, the Hessian H of V is modulated a slow time scale with f ((vo))) proximity of the fixed points Vaa. We call A the effective Hessian matrix.
15

---

Published as a conference paper at ICLR 2025
$1.4 EFFECTIVE HESSIAN
Stability of the fixed points is determined by the eigenvalues of the Jacobian matrix J. The blocks Jea and Jan. are such that (see also (Leleu et al. 2019)):
Jew bl,
(S17)
with I the identity matrix and -5aa vs(va) and Jee 0. Consequently, the determinant 1/ (/((a)
P(A) -| A characterizing the eigenvalues of J has the following property:
Jax A1 Jer
Jue -X = -ax - ΛTIA JacJez = =1az NIA - bT.
P(A)
(S18)
we note μ4 the eigenvalues of Jaw the characteristic polynomial can be rewritten:
P(A) - A)^ bt X² μd
(S19)
The eigenvalues X can come in pairs X and x depending on the sign of A, = 11 + and are expressed as follows by solving P(A) = 0 (Aj E .N):
x= = (ug + A 2 if i > 0,
(S20)
ス = (m) / A 0.
(S21)
We define the pre-conditioned Hessian matrix A as Al ((()) = D[ pl o.
H (Vao,) such that
ozū (a.va)
DicaM(Vao:)H DI
H H(Voo),
(S22)
av D2.
(aiva)
v
D
HH(vao)
(S23)
Gz. Va(Va)
(oiva)
1 H(Vao),
(S24)
v
The last step uses the fact that
(z)) (x)
is odd function.
Note that the factors
Vi, multiplying the Hessian H positive at local minima satisfying
av Der (d.va)'
OV Dz, (a. ViVaa, O, Vi. In this case, the eigenvalues of D (oiva) H[ are the same as D 1/2H(Voo;) ( "ل 1/2 By Sylvester's s law of matrix inertia, this implies that
deg
the eigenvalues of Á and H have the same signs. matrix theory, the term "inertia" refers to the number of positive, negative, and zero eigenvalues of real symmetric matrix.
av Lemma 2: The signs of the eigenvalues of Å and H are the same when a; dr, (a.va) > O, Vi.
Theorem The effect error variables at proximity of fixed points to replace the Hessian H of V by an effective Hessian / with H (Voor Dí av dap (a:va) [H( The inertia of the Hessian and
0V (aiva) Vi. effective Hessian are the same for local minima of V verifying the condition a; Bz,
16

---

Published as a conference paper at ICLR 2025
$1.5 EIGENVALUES H AND STABILITY OF THE FIXED POINTS
The eigenvalues of the Jacobian Jez- noted (4j. can be expressed using the eigen spectrum of the pre-conditioned Hessian H, noted 7/. as follows using eq S7
+ a Va (/((a) 74
Pj -af(Vao))
(S25)
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. S20 [220 and S2T
(Rel( + /2z + 46) = 0 pf 4b 0,
(S26)
(Reljj + iv 11 4b)) = 0 if n3 4b < 0.
(S27)
Noting 7 the eigenvalue with maximal part Re Re/n,), eqs.
S26 and S26 imply that the fixed point Va
becomes stable under the condition (see eq. S21
0,
(S28)
@vV i... - af((aa,,) eig(Voo,) 5 0.
(S29)
Assuming 0x? 0, we have (see eqs. S3 and S12
(x())) OV onbg(a))) <g((;;)
(S30)
(2,(,,))) dz, oV de,
a (1% (xx))
(S31)
Moreover, have (see eq. SI11:
4% (m:)
f(x,)
(S32)
Oxi
= 2
(S33)
(x:())
Combining eqs. S31 and S33 eq. S29
can be written as follows:
7(Vao) = 2 1/ (V(
(S34)
1/ (va)
Va(a
Theorem For small values the target amplitude a such that ((Vao) > 2 the cigenvalue of Ai ((ar) with (VTor) = D
Va(Va) where ((vao)
(o)o H ((vao,)) the system described by
eqs. $2 and S3 does not have any stable fixed points assuming
@VV 0.
8x
17

---

Published as a conference paper at ICLR 2025
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(Φ) = = Eij Wiyzyaa the calculations are simplified because av (a./a) va 0V dz. (o:) and (a (aV and eq. S34 can be written:
(() 2 11 ( a
(S35)
(
where ÿ is independent of Va.
18

---

Published as a conference paper at ICLR 2025
a
1.0
P=P(aa)
X
7a=088
0.5
0.60.6
0.0
0.5
1.0
1.5
2.0
3.5
4.0
F(ya)
b
0
0
20
40
60
80
100
20
40
60
80
100
d
2
20
40
60
80
100
0
20
60
80
100
Figure S1: (a) Bifurcation diagram of the system of eqs. S2 and S3 in the case of the V(x) Ei Wey0;0 with = 8 and Way is randomly chosen to be equal to or -1 for j > í. Wij = wij and Wii = O. The red curve corresponds to the set {vā, F( ((va) with > 0 which denote the condition under which the real part of the eigenvalues μ of the Jacobian matrix JaN is zero. The and X markers correspond to the maximum real part 7 of the cigenvalues of DI pr (..vāj ( (Vaoo)
SE;
for the x-axis and the value of 0.6 and a 0.8, respectively. (b) and (c) Time series of soft spins t; vs. time in the case 0.8 and a = 0.6, respectively. (d) and (e) The same as (b) and (c) for the error variables e;.
In Fig. SI the stability of the dynamics described in eqs. S2 and S3 when Ising Hamiltonian with N 8 is shown by the bifurcation diagram in the space F
va). The + and x markers (a.va) ((vao,) for the x-axis
correspond to the position real part of the eigenvalue
and the value of 0.6 and defined by [Va.F((a) ((va) with a > 0 (shown red). this regime, the system settles to a stable fixed point as shown in Fig. [1 and (d). For 0.6, the point {5, v is below the red curve. In this regime, the system does not possess any fixed points and the limit behavior a periodic cycle.
0.8, respectively. For 0.8, the point {7, va is above the set
19

---

Published as a conference paper at ICLR 2025
S2 PSEUDO-CODE OF MHCACM
The Pseudo-code of Cetropolis-Hastinng adjusted chaotic amplitude control with momentum (MH- CACm) is shown in this section.
Algorithm 1 Hattsiso adjusted chaotic amplitude control with momentum
function DETERMINISTIC PATH())
Iterate deterministic steps CACm
2:
u
< 0 3:
4:
"" <-0
0
5:
6:
y t-a
Initialize
7:
for I f- to n do
8:
u" f u'
9:
u' u on'' (V) +7(1' -0''=""))
Fuler steps
et - olyoy a) oe
12:
et <0
Normalize
13:
f Ite Fa
14: end for
15: y Ltt 20 (aa)
16: return
17:
end function
18: function PROBABILISTIC JUMP(y)
Generate random jump
19: T= with probability 2=1
2 2 - otherwise
20: return
21: end function
22: function METROOPOLISHAASINNEE STEP(T, o, R.
Metropolis- Hasttingss- step
23:
A(T/o) <- min(1, P(a)Qt//a)) F.,,
P(riQ(a/T))
Di Function F defined in (9)
24: B TRUE with probability A(T/o), B FALSE otherwise
25:
if B then
Accept new configuration
27:
0 T
28:
end
29:
return T.
30:
end function
procedure SAMPLE (V)
32:
Y < 440,, 1)
D. Random initialization
33:
a + PROBABILISTIC JUMP (v)
> Initial configuration
34:
for t. K de
Generate K samples
PROBABILISTIC JUMP (y)
New candidate configuration
36:
DETERMINISTIC PATH (T)
37:
0. < METROPOLIS- HITRROOOOISSSSSTINGSS HASTIGGS STEP (T,
38:
end for end procedure
20

---

Published as a conference paper at ICLR 2025
S3 BOLTZMANN SAMPLING
S3.1 VERIFICATION BOLTZMANN SAMPLING AT EQUILIBRICM
To test our Markov chain approach, we sample from the Boltzmann distribution of the Ising Hamilto- nian (o = Eij Wijo,oj. We use the Wishart planted ensemble (WPE) (Hamze et al. 2020 of size N 18 spins to generate the Ising coupling matrix. Using the WPE to generate instances the Ising Hamiltonian allows us to tune the complexity of the energy landscape and also know the ground state a-priori (Hamze et al. 2020). 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. 2020). 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 3 = 0.07 by calculating the symmetric KL divergence (see Appendix section (3.2). Note that CACm (without the sample from the target Boltzmann distribution (see Fig. S2 a). In contrast, MHCACm exhibits exactly the same distribution of energy as the Boltzmann distribution thermal equilibrium. This is accomplished while retaining the global exploration properties of CACm, making it potent hybrid heuristic for identifying lower energy states of non-convex energy functions.
step) does
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(ar = Eij Wyyra for
Wishart planted instance (Hamze al., 2020) with a size of N 18 spins an inverse temperature of 3.
CACmn
1.4
CACm MIICACm
MECACCm (proposed)
Gibbs distribution
1.2
KL divergence
inv. compenature
1.0
pdD
10
0.06
CMOs, 20260 MICAAC D/Y
0.01
MICACE
0.2
04
0.02
0.03
0.04
0.05
21
X101
x- Punfff
effective inv. rempertture Par
Figure S2: (a) Exact Boltzmann distribution (green dotted line), sampled distribution using CACm (without MH step, red line), and MHCACm (without MH step, see blue line). Gaussian filtering of the energy distribution is applied for the sake of readability. , = 0.07. (b) KI. divergence vs. ratio of inverse temperatures λ =8- the KL divergence is minimized vs. B. The states used for the sampled distribution are boolean configurations accepted by the MH criterion.
Par (c) Effective inverse temperature & temperature at which
S2 is shown that the KL divergence is minimized for a ratio of inverse temperatures Note that Ber is the effective temperature defined in eq. of the main manuscript. The effective inverse temperature Ber at which the divergence is minimized temperature in noted & Figure S2 (c) shows that 3 is linear function of 3 within the range of parameters considered. Moreover, the addition of the MH step to MHCACm increases the similarity between 8 and B compared to that of CACm.
In Figs.
Bur 8
S3.2 SIMULATIONS FOR SMALL N
In Fig. S2 the distance the Boltzmann distribution P(o) (V))) (Z= E. Al(() ) is calculated
using the symmetric KL divergence. The symmetric KL divergence D (P Q) between the distributions and Q defined as follows:
21

---

Published as a conference paper at ICLR 2025
1 DskL 1) DKL (P (1 Dx Q PJ),
(S36)
where
P(i) Dxi (p Σ P()) log Q(1
(S37)
In the case of Fig. S2 exact Boltzmann distribution of energy levels obtained by brute force calculations and the estimated adjusted chaotic
the convergence of the MCMC is measured by the KI. divergence between the
distribution of energy levels obtained by simulation of the Metropolis-Hastinn discrete space sampling. The former computation is possible because the problems size N = 18 sufficiently small for the brute force search method. The instance used for Fig. 2 randomly generated instance from the Wishart planted set with hardness parameter owwPF 0.8 (Perera et al. 2020).
S3.3 SIMULATIONS FOR LARGE N
In order test the convergence of the MCMC in the case of larger N. utilize heuristic metric. Conventional Gibbs sampling algorithm is run for a sufficiently large number of steps for the chain converge. In the case of N 60 and owpe 0.8, numerical simulations indicate that the convergence occurs within about 1000 time step of the Gibbs sampling algorithm. To observe the convergence, limit the characterization of the Boltzmann distribution to its first two moments: the mean energy < H > with < H >= E. H(ar) P(o) and standard deviation of the energy V (H H > Figure S3 shows the convergence of the mean and standard deviation of the energy under Gibbs sampling.
x10
1.2 X10
1.0
217>
-2
0.6
-3.
0.4
4-
0.2
-5
200
400
800
1000
0.0
600
1000
Monte Carlo steps
Monte Carlo steps
Figure S3: Convergence of the mean and standard deviation of energy distribution of a N = 60 Wishart planted instance with O'we 0.8 under Gibbs sampling. Approximate convergence is reached after 1000 Monte Carlo steps. $ 2.5.
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 then calculated as follows:
D Hehaos >
HGibbs
>
+ (Heuus
H >>
V< (Hobb.
H >>>>>
(S38)
22

---

Published as a conference paper at ICLR 2025
where the subseript denote the cases of Gibbs and MH chaotic sampling, respectively.
In Fig. $4 is shown distance D vs. the number of steps T in the algorithm. One step corresponds to 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 the steady-star distribution. The of auxiliary variables speed up the convergence compared to the case without.
107
Pripeesi g
Propored
Gilis compling
105
Distance and
10*
10a
10a
10
T
Figure S4: Distance D between the mean and variance the distribution obtained by the proposed Metropolis-Hastinngs adjusted chaotic discrete space sampling and the obtained by iterating the Gibbs sampling algorithm (see eq. S36) vs. the number of steps in the algorithm. One step
corresponds to tentative spin flip and one Euler step in the Gibbs sampling and chaotic sampling case, respectively.
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. 2020). Numerical simulation the convergence 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 the use of the replica exchange Monte Carlo method in both cases of Gibbs sampling and proposed chaotic sampling.
23

---

Published as a conference paper at ICLR 2025
S4 AUTOTUNING USING DYNAMIC ANISOTROPIC SMOOTHING
The parameters of the solvers are described in Table ST We have run the autotuning algorithm called dynamic anisotropic smoothing (DAS) (Reifenstein al. 2024) for each problem setting (problem size N, bias b). The MH sampling rate is set to JMM 0.1.
Table S1: Parameters of the algorithms tuned using DAS
Algorithm
number of parameters
parameter list
SA
3
Bii Bän T
AIM
a, 7, T
CAC
5
a, T
CACm
B, a, E, a, 7/4 Ni O, &, a, 71
HMCACM
An example of the the evolution of parameters during the DAS autotuning shown in Fig. S5
beta
kappa
100
lamb
xi
gamma
parameters
10-1
0
50
100
150
200
250
300
steps
Figure S5: Evolution of parameters during the DAS autotuning for MHCACm. N 100, b 12, T' 1000. Here "lamb" denotes a and "kappa" is A
3
24

---

Published as a conference paper at ICLR 2025
S5 CONTRUCTION OF DEGENERATE WISHART PLANTED ENSEMBLE
The Wishart planted ensemble instances are generated by first choosing a planted ground state g E [-1,1)% However, our case we consider of planted ground states which are represented as columns an N by matrix G E (-1,11)XD, Let & be a matrix with the same columns space as G whose columns are orthonormal. Then we generate an N by M matrix W of i.i.d. Gaussian random variables. This matrix is then transformed so that the columns orthogonal to all of the ground states form the matrix given follows: W
GG W.. The Ising coupling matrix constructed as j WWT with the diagonal elements being set to zero. We can see that every column of G is ground state since x Jæ 2 0 and G JG 0.
Note that in the Wishart planted ensemble, the parameter awpE N where N is the num- ber rows and M is the number of columns, quantifies the aspect ratio of the Wishart matrix used to define the interaction strengths the Ising model. This parameter plays a central role in determining the ruggedness of the energy landscape.
To add bias the energy landscape, 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 = GG GAW where A is symmetric matrix and W is the random matrix defined above. In this work we consider A to be defined A = I ++ $ 11 where is the "bias" away from the ferromagnetic solution. The point this that the residual energy (energy away from the ground state) will be amplified for states that are close the ferromagnetic solution. Then, when choosing the planted ground state can choose one which is close to the ferromagnetic solution and one that is far creating bias in the search towards ground state with a hamming distance of 3 from the ferromagnetic solution and one ground state chosen completely randomly from { {-1, 1] A Gauge transform is applied obfuscate the ground states at the end of this process.
somehow create the instance such that
ground state another. For the numerical results in this work consider one
S6 COMPARISON TO RECENTLY PROPOSED SAMPLING ALGORITHM
In this section, we compare the proposed algorithm to the sampling method described in (Sun et al. 2023). We have implemented the algorithm of (Sun et al. 2023), that we denote by the acronym ISCO, and compare it against MHCACm in table $2. Success probability is higher using MHCACm.
Po (N=60)
iSCO
MHCACm
b=0 (no bias)
0.0003
0.0008
b= 12 (bias)
0.51
0.15
Table S2: Probability to find any ground-state from unbiased (b 0) and biased (b = 12) Wishart planted instances using iSCO and MHCACm after 500 time steps. N
S7 COMPUTING ENVIRONMENT
All simulations are done in CPU using a Intel Core i9- 19-11H 8 cores, 2.6 GHz and GPU using 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.
25


Reviews:
[{"review_id": "ddwyrKo0YU", "reviewer": "Reviewer_ymkA", "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.", "questions": "NA", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The empirical results include multiple baselines and comparisons in terms of different metrics.\n- The visualization in Fig.1 clearly shows the main algorithmic idea.", "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.\n- 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? \n- 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.\n- The runtime comparison only considers simulated annealing. It will be better to include other baselines as well.\n- 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].\n\n[1] Revisiting sampling for combinatorial optimization, ICML 2023", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8Wml87G7nb", "reviewer": "Reviewer_wvMs", "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.", "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.\n\n(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.\n\n(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. \n\n(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.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 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.\n\n(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.", "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. \n\n(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. \n\nIn 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "pyIl0mrOG0", "reviewer": "Reviewer_g5qM", "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.", "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?\n* 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.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The paper writing is nice and structured, with a comprehensive literature review and detailed problem set-up.\n* The paper looks technically sound with solid mathematical proof.\n* The proposed algorithm is evaluated on multiple tasks and compared with various other benchmark methods, showing competitive performance.", "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.\n* The connection between ground-state sampling and deep learning optimization/generalization mentioned in the paper is interesting, but the discussion is very limited. \n* For numerical experiments, the definition of TTS is hard to comprehend. Does smaller TTS indicate better algorithmic performance?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Fg9wQehEn2", "reviewer": "Reviewer_zbqu", "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.", "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. \n- Can you elaborate more on the \"dual-primal Lagrangian approach\" and its difference from CAC?\n- Some minor typos: e.g., line 122 \"in order (to) benchmark this algorithm's ability\".", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 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. \n- 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. \n- MHCACm is well-suited for large-scale deployment on GPU because its computational bottleneck is matrix-vector multiplication.", "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.\n- 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.).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8udiuizOT8", "reviewer": "Reviewer_SwgT", "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.", "questions": "See Weaknesses.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 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.\n2. The hybrid algorithm shows potential for improving ground-state sampling performance in combinatorial optimization\n3. The authors address parallelization for efficient computation on GPUs\n4. 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.", "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.\n\n\n**Major concerns**\n1. 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.\n2. The paper lacks theoretical guarantees regarding the convergence or performance of the proposed algorithm.\n3. 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}$.\n\n**Other suggestions**\n1. The paper lacks clear definitions for essential variables (e.g., $\\boldsymbol{u}$ and $\\boldsymbol{e}$ in Equation (1)).\n2. The time variable $t$ is used ambiguously. It denotes continuous evolution in Equation (1) but has discrete updates in Equations (4)-(5).\n3. 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.\n4. 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)