Efficient Combinatorial Optimization via Heat Diffusion
Authors: Hengyuan Ma
Wenlian Lu
Jianfeng Feng
Venue: NeurIPS 2024 poster (2024)
Keywords: Combinatorial optimization
Heat equation
Primary area: optimization
TL;DR: We propose a new framework for combinatorial optimization based the insight that heat diffusion can propagate information from distant area to the solver for efficient navigation.

Abstract: Combinatorial optimization problems are widespread but inherently challenging due to their discrete nature. The primary limitation of existing methods is that they can only access a small fraction of the solution space at each iteration, resulting in limited efficiency for searching the global optimal. To overcome this challenge, diverging from conventional efforts of expanding the solver's search scope, we focus on enabling information to actively propagate to the solver through heat diffusion. By transforming the target function while preserving its optima, heat diffusion facilitates information flow from distant regions to the solver, providing more efficient navigation. Utilizing heat diffusion, we propose a framework for solving general combinatorial optimization problems. The proposed methodology demonstrates superior performance across a range of the most challenging and widely encountered combinatorial optimizations. Echoing recent advancements in harnessing thermodynamics for generative artificial intelligence, our study further reveals its significant potential in advancing combinatorial optimization.

Efficient Combinatorial Optimization via Heat Diffusion
Hengyuan Ma
Institute of Science and Technology for Brain-inspired Intelligence
Fudan University
Shanghai, China 200433
hangyuanma210m fudan edu. cn
Wenlian
Institute of Science and Technology for Brain-inspired Intelligence Fudan University
Shanghai, China 200433 wenl i an@fudan edu. cn
Jianfeng Feng
Institute of Science and Technology for Brain-inspired Intelligence
Fudan University
Shanghai, China 200433 j jian eng640gmail com
Abstract
Combinatorial optimization problems are widespread but inherently challenging due their discrete nature. The primary limitation of existing methods is that they can only access small fraction of the solution space at each iteration, resulting in limited efficiency for searching the global optimal. To overcome this challenge, diverging from conventional efforts of expanding the solver's s search scope, we focus on enabling information actively propagate the solver through heat diffusion. By transforming the target function while preserving its optima, heat diffusion facilitates information flow from distant regions to the solver, providing more efficient navigation. Utilizing heat diffusion, we propose framework for solving general combinatorial optimization problems. The proposed methodology demonstrates superior performance across range of the most challenging and widely encountered combinatorial optimizations. Echoing recent advancements harnessing thermodynanics for generative artificial intelligence, reveals its significant potential advancing combinatorial optimization. The codebase of our study is available in battps://githee com/ AwakerMhy/Hee
study further
Introduction
Combinatorial optimization problems are prevalent in various applications, encompassing circuit design [ machine learning (2), computer vision [3], molecular dynamics simulation [4), traffic flow optimization [5]. and financial risk analysis
(6]. This widespread application creates a significant demand for accelerated solutions to these problems. Alongside classical algorithms,, which encompass both exact solvers and metaheuristics (7), recent years have seen remarkable advancements in addressing combinatorial optimization. These include quantum adiabatic approaches [8 9 I0]
simulated bifurcation [ 12 13), coherent Ising machine [4 [5]. high-order Ising machine [6]
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

and deep learning techniques [17 [8]. However, due to the exponential growth of the solution number, finding the optima within a limited computational budget remains a daunting challenge.
Our primary focus is on iterative approximation solvers, which constitute significant class of combinatorial optimization methods. An iterative approximation solvers typically begin with an initial solution and iteratively improve it by finding better solutions within the neighborhood of the current solution, known as the search scope or more vividly, receptive field. However, due to combinatorial explosion, as the scope of the receptive field increases, the number of solutions to be assessed grows exponentially, making thorough evaluation of all these solutions expensive. As result, current approaches are limited to narrow receptive field, rendering them blind to distant regions the solution space and heightening the risk of getting trapped in local minimas areas with bumpy landscapes. Although methods like large neighborhood search [19], variable neighborhood search 20 and path auxiliary sampling [2]] are designed to broaden the search scope, they can only gather a modest increment of information from the expanded search scope. Consequently, the current solvers' receptive field remains significantly constrained, impeding their efficiency.
original problem
problem under heat diffusion
time
high
Where can find the key?
flow over here.
know where
to go next!
temperature
key
low
cooperative optimization
Figure 1: The heat diffusion optimization (HeO) framework. The efficiency of searching a key in dark room is significantly improved by employing navigation that utilizes heat emission from the key. In our framework, heat diffusion transforms the target function of a combinatorial optimization problem into different versions while preserving the location of the optima. Therefore, the gradient information of these transformed functions cooperatively help to optimize the original target function.
In this study, we approach the prevalent limitation stated above from unique perspective. Instead of expanding the solver's receptive field to acquire more information from the solution space, we concentrate on propagating information from distant areas of the solution space the solver via heat diffusion (22). To illustrate, imagine a solver searching for the optima in the solution space akin to a person searching for key dark room, as depicted Fig. Without light, the person is compelled to rely solely on touching his surrounding space. The tactile provides only localized information, leading inefficient navigation. This mirrors the current situation in combinatorial optimization, wherein the receptive field is predominantly confined to local information. However, if the key were to emit heat, its radiating warmth would be perceptible from a distance, acting as directional beacon. This would significantly enhance navigational efficiency for finding the key.
Motivated the above metaphor, we propose simple but efficient framework utilizing heat diffusion to solve various combinatorial optimization problems. Heat diffusion transforms the target function into different versions, within which the information from distant regions actively flow toward the solver. Crucially, the backward uniqueness of the heat equation [23] guarantces that the original problem's optima are unchanged under these transformations. Therefore, information of target functions under heat diffusion transformations of different degrees can be cooperatively employed for optimize the original problem (Fig. [].. Empirically, our framework demonstrates superior performance compared to advanced algorithms across a diverse range of combinatorial optimization instances, spanning quadratic to polynomial, binary to ternary, unconstrained to constrained, and discrete to mixed-variable scenarios. Mirroring the recent breakthroughs in generative artificial intelligence through diffusion processes www phenomenon, in enhancing combinatorial optimization.
24 our research further reveals the potential of heat
diffusion, related thermodynamic

---

2 Failure of gradient-based combinatorial optimization
We will first formulate the general combinatorial optimization problem and then reframe gradient descent-base bascent-based manner. This reformulation allows us to utilize heat diffusion later. Various combinatorial optimization problems can be naturally formalized as a pseudo- Boolean optimization (PBO) problem [25], in which we aim to find the minima of a real-value target function f E I + R subjecting to binary constraints
min [(s) SE ==(-1,1"
(I)
where is binary configuration (bits), and I(.) is the target function. Through the transformation
+ 1)/2, our definition aligns with that [26] where elements of s take 0 or Given the advanced are interested in converting the discrete optimization problem into a differentiable one, thereby enabling gradient descent. To achieve this purpose, we encode the bits St, i I,..., independent Bernoulli variables p(s; = 0.5 = (0; 0.5)
development of gradient-based algorithms,
with 0g E [0, 1J. In this way, we convert the original combinatorial optimization problem into
min (0), OET
(2)
where I := [0, 1]", and h(0) = Eptsto) The minima 0 of Eq. is 0 (.5(sgn(s*) + 1), given that s* is minima of the original problem Eq. (I.. Here, sgn(.) is the element-wise sign function. Now Eq. (2) can solved through gradient descent starting from a given initial Oo
0 0: Ve h(0.),
1=
where is the learning rate and T is the iteration number. Unfortunately, this yields probability distribution p(s(0) over the configuration space [-1, instead of a deterministic binary config- uration desired. Although we can manually binarize 0 through ((0) == sgn(0 0.5) to get the binary configuration which maximizes probability p(s[0), the outcome (() may be much higher than h()), resulting in significant performance degradation [27]. This suggests that good gradient-basse baased optimizer should efficiently diminish the uncertainty in the output distribution p(s(0). which can be measured by its total variance
v = Σ 0,(1-0:) (c)- 02).
2.1 Monte Carlo gradient estimation
Conventionally, solve the problem Eq. by approximating the gradient of h(0) via Monte Carlo gradient estimation (MCGE) (28) (Alg. Appendix). which we estimate the gradient in
Eq.
Veh(O) (p(s))))((s)) vo (og)())()) ~ M Σ ((s(m)) Vo log (m)
00)
where s(m) m..idd p)s m = 1,.... M. However, turns out that MCGE performs poorly even equipped with momentum, compared to existing solvers such simulated annealing and Hopfield neural network, as shown in Fig. 2 We interpret this result as follows. Although MCGE turns the combinatorial optimization problem into differentiable one, does not reduce any inherent complexity of the original problem, which may contains a convoluted landscape. As gradient only provides local information, MCGE is susceptible be trapped in local minimas.
3 Heat diffusion optimization
The inferior performance of the MCGE is attributed to the narrow receptive field of the gradient descent. To overcome this drawback, we manage to provide more efficient navigation to the solver by employing heat diffusion 29). which propagates information from distant region to the solver. Intuitively, consider the parameter space as thermodynamic system, where each parameter referred to location and is associated with an initial temperature value -h(0), shown in Fig,
3

---

Then the optimization procedure can be described as the process that the solver is walking around the parameter space to find the location 0r with the highest temperature (or equivalently. the global minima of I())). As time progresses, heat flows obeying the Newton's ' law of cooling, leading to an evolution of the temperature distribution across time. The heat at 0" flows towards surrounding areas, ultimately reaching the solver's location. This provides valuable information for the solver, as it can trace the direction of heat flow to locate the 0".
3.1 Heat diffusion on the parameter space
Now we introduce heat diffusion for combinatorial optimization. We extent the parameter space of 0 from [0, 1 to I with R R (-x0, +oo}. To keep the probabilistic distribution p(s[0)
meaningful for O & [0, 1]", we now redefine p(s; = =1 0) = clamp(0.5 + (0; - 0.5), 0, 1), where
the clamp function calculated as caamp(æ; 0, 1) = max (0, (i. 1)). Denote the temperature at
location & and time r as u(T, 0), which is the solution to the following unbounded heat equation [29]
(,,(TT,)) 0) = Dou(T, 0), u(TT O) = 1()),
>> 0, 7= 0,
i E R R
(6)
where A is the Laplacian operator: Ag(x) Ei= 8z (x). For 0 E RR"/R", we define M(T, 0) limg,,-0 u(T, 0,), where (0..) is a sequence T8" converged to 0. Heat equation in the combinatorial optimization exhibits two beneficial characteristics. Firstly, the propagation speed of heat is infinite [22), implying that the information can reach the solver instantancously. Secondly, the location of the global minima does not change across time т, demonstrated in the following theorem.
Theorem 1. For any the function u(T, 0) and h(0) has the same global minima in I arg mingge u(T 0) arg mingge=.h(())
Consequentially, we can generalize the gradient descent approach Eq. (3) by substituting the function h(0)) with u(Tr, 0) for different > 0 at different iteration step t in Eq. (3. as follows
0.+1 = 0. - 7 Vo u(TT, Ot),
where the subscript ' in T, means that Tr can vary across different steps. In this way, the solver can receive the gradient information about distant region of the landscape that is propagated by the heat diffusion, resulting in more efficient navigation. However, Eq. (8) will converge to
0. +00 project the 0, back to I after gradient descent each iteration
since O, are unbounded. To make the procedure Bq. { practicable,
00, $f
0.+1 = Projz (0: 7 Vo Ot)).
(9)
so that Og E I always holds, where we define the projection as Projz(x); = min(1, (1, max(0, x;)) for I,..., and for E E". Eq. (9) reasonable update rule for finding the minimum 0* within I for the following two reasons: (1) The projection of the Ó is the minimum of h(() in I, i.e., Projz (0) 0": (2) Due to the property of the projection, if the solver moves towards 0, it also gets closer to 0" More importantly, since the coordinates
are all infinite, the convergent point of Eq. (9) must be one of the vertices of I, i.e., (0, 1 1/n, This suggests that Eq. 9 tends to give an output that diminishes the uncertainty (()) (Eq.
3.2 Solving the heat equation
To develop an algorithm for solving the combinatorial optimization problem from Eq. 99) we must solve the heat equation Eq. (6). which seems a significant challenge when the dimension high. Fortunately, the solution has a closed form if the target function [s( can be written as a multi-linear polynomial
(s) ao + Σ 81, +
Σ 113
i2
aKJj..ig
Sige
(10)
condition met by a wide range of combinatorial optimization problems [66]

---

Theorem 2. Supposed that f(s) is a multilinear polynomial of's, then the solution to Eq. 6
WTT 0) = pex (f
)
x x
(1))
where erf( 2 /x So e dt the error function.
For more general cases other then Eq. (00) (such as Eq. (Eq.
13
we still use the approximation
3.3 Proposed algorithm
Based Eq. 9 and Thm. proposed heat diffusion optimization (HeO), saadletttbbsed algorithm for combinatorial optimization, as illustrated in Alg. [ where we estimate Eq. (II) with one sample X, Unif(0, 1J", and we denote VT, ot for short. Our HeO can be equipped with momentum, which is shown in Alg. i Appendix. In contrast to those methods designed for solving special class of PBO (Eq. (T) such quadratic unconstrained binary optimization (QUBO), our HeO can directly solve PBO problems with general form. Although PBO can be represented
as QUBO [30. this necessitates the introduction of auxiliary variables, which may consequently increase the problem size and leading to additional computational overhead [T]] Compared to other algorithms, our HeO has relatively low complexity. The most computationally intensive operation at each step is gradient calculation, which can be explicitly expressed or efficiently computed with tools like PyTorch ss autograd, and can be accelerated using GPUs. As shown in Fig. of the in Appendix, the time cost per iteration of our methods increases linearly with the problem dimension, with a small constant coefficient. Therefore, our HeO efficient even high-dimensional cases.
Algorithm 1 Heat diffusion optimization (HeO)
Input: target function f-), step size Y, a schedule { iteration number T initialize elements of 0, as 0.5 for t to T 1 do sample x ~ Uni 1j
g: < Vol 0 xx=x))
0--1 - Projz ( - >ge)
end for
Output: binary configuration ST sgn(Or 0.5)
One counter-intuitive thing is that to minimize the target function h(0), the HeO actually are minimizing different functions u(Tr, 0) at different step t. We interpret this by providing an upper bound for the target optimization loss h(0) h(0)).
Theorem 3. Denote f max, fs)). Given T2 O and 0, there exists E (0, such that
1(0) h(00))))))0 < [(] (') u(72. 0*) + u(T, 0) u(T, 0-) dr) 1/2 + (12)
T
Accordingly, minimizing u(r, 0) for each cooperatively aids in minimizing the original target function h(0). Thus, we refer to HeO as a cooperative optimization paradigm, as illustrated in Fig.
Experiments
We apply our HeO to variety of NP-hard combinatorial optimization problems demonstrate its broad applicability. Unless explicitly stated otherwise, we employ the 7; schedule as / =: a; vē(1 - t/TT) for HeO throughout this work. The sensitivity of other parameter settings including the step size y and iterations are shown in Fig- S2 This choice is motivated by the idea that the reversed direction of heat flow guides the solver towards the original of its source, i.e., the global minima. Noticed that this choice not theoretically necessary, as elaborated in Discussion.
Toy example. We consider the following target function
f (s) = a. sissiggoi +ai)
(13)
5

---

SA (S steps)
energy
SA (5 10" steps)
22000
(5 10 steps)
SA (5 10 steps)
SA (5 steps)
MCGE
-2000
ICGE//mmmenumm)
HNN
uncertainty
2000
HeO
1000
10F
10a
10?
107
10*
10a
10f
107
10a
time step
Figure 2:
Performance of HeO (Alg. Appendix), Monte Carlo gradient estimation (MCGE), Hopfield neural network (HNN) and simulated annealing (SA) on minimizing the output of neural network (Eg. (3.). Top panel: the target function. Bottom panel: the uncertainty V(0) (Eq. (4)).
grapth
77/
6.002%
594
3.70%
45%
2.00114
osa examples of
1/66
0.00%
1.10%
0.611
2%:
3.54%
1.20/4
5.00%
0.00%
0.22%
116
5.36%
3.20%
70
ostrunber: CA Sumber
40
75
27
wormppe of
50
36-
30
20-
28-
25
15-
10
14
160 100
B5R
358 R60
ESA dsa
He0 10A
Figure 3: a, Illustration the max-cut problem. Performance of HeO (Alg. and representative iterative approximation methods including LQA [10), aSB [12. bSB [13. dSB [3, CIM (35] and SIM-CIM [15] on max-cut problems from the Biq Mac Library (36). Top panel: average relative loss for each algorithm over all problems. Bottom panel: the count of instances where each algorithm ended up with one of the bottom-2 worst results among the 7 algorithms.
where sigmoid(z) 1 and the elements of the network parameters a1 R" ag Rm W Rm.. uniformly sampled from [-1, 1] and fixed during optimizing. According to the universal approximation theory [32). f(s) can approximate any continuous function with sufficiently large m. thereby representing general target function. We compare the performance of our HcO with momentum (Alg. 3) against several representative methods: the conventional gradient- based solver MCGE 28 (with without momentum), the simulated annealing [33,, and the Hopfield neural network 34 As shown in Fig. 2 our HeO demonstrates exceptional superiority over all other methods, and efficiently reduces uncertainty compared MCGE.
Quadratic unconstrained binary optimization (QUBO). QUBO is the combinatorial optimization problem with quadratic target function (J E paan symmetric matrix with diagonals)
f(s) = ST Js.
(14)
This corresponds to the case where K = in Eq. 10).. A well-known class of QUBO max-cut problem (27]. in which we divide the vertices of a graph into two distinct subsets and aim to maximize
the number of edges between them. Its target function is expressed as Eq. (14). where determined by the adjacency matrix of the graph.
We compare our HeO with representative iterative approximation methods especially developed for solving QUBO including LQA [10), aSB [2), bSB [13]. dSB [33) CIM [35], and SIM-CIM [15]

---

Congundiive normal form (GNF)
1.00
HeO
bi thy li bul the yid ba bil
DE [0.3)" 33"
5-1-b
eat-sted conetrinnt
Vi Set sftee probtbi byy
Circuil
D
OR
D
AND
0.98
sallidying or
number varisbles
number of variables
Figure 4: Illustration of the boolean 3-satisfiability (3-SAT) problem. b, Performance of HeO (Alg. [4, Appendix), 2-order and 3-order oscillation Ising machine (OIM) [16] on 3-SAT problems with various number of variables from the SATLIB [37. We report the mean percent of constraints satisfied (left) and probability of satisfying all claims (right) for each algorithm.
max-cut problems in the Biq Mac Library 36 We report the relative loss averaged over all instances and the count of the instances where each algorithm gives the bottom-2 worse result among the 7 algorithms. As shown in Fig. 3 our HcO is superior to other methods in terms of two metrics.
class of combinatorial optimization problems, in which higher-order interactions between bits appears in the target function. Existing methods for solving PUBO fall into two categories: the first approach involves transforming PUBO into QUBO by adding auxiliary variables through a quadratization process, and then solving as QUBO problem (88. and the one directly solves PUBO [16]. Quadratization may dramatically increases the dimension of the problem, hence brings heavier computational overhead, while our HcO can be directly used for solving PUBO. well-known class of PUBO is the Boolcan 3- satisfiability (3- (ATS problem (27). which involves determining the satisfiability of Boolean formula over n Boolean variables bi..... b, where 6; E {0, 1}. The Boolean formula is structured in Conjunctive Normal Form (CNF) consisting of H conjunction (A) of clauses, and each clause h is disjunction (V) of exactly three literals (either Boolean variable or its negation). An algorithm of 3-SAT aims to find the Boolean variables that makes many clauses satisfied.
Polynomial unconstrained binary optimization (PUBO). PUBO
To apply our HeO to the 3-SAT, we encode each Boolean variable by as 8,, which assigned with value 1 if b, 1, otherwise = -1. For a literal, we define value CR. which is -1 if the literal the negation of the corresponding Boolean variable, otherwise variables that satisfies many as clauses is equivalent minimize the target function
Then finding the Boolean
f(s)
ΣI
1
Ch, 2
(15)
A=
This corresponds to the case where K = 3 in Eq. 1T We compared our HeO (Alg 4 Appendix) with the second-order oscillator Ising machines (OIM) solver that using quadratization and the state-of-art 3-order OIM proposed in [16] on 3-SAT instance in SATLIB As shown in Fig. 4 our HcO is superior to other methods in attaining higher quality solutions and finding more the complete satisfiable solution (solutions that satisfying all clauses). Notably, for the cases of 175-250 variables, our HeO able to find more complete satisfiable solutions, compared to the 3-order OIM, while the 2-order OIM fails find any complete solutions [16).
Ternary optimization. Neural networks excel learning and modeling complex functions, but they also bring about considerable computational demands due to their vast number of parameters. promising strategy to mitigate this issue is quantization, which converts network parameters into discrete values [3)]. However, directly training networks with discrete parameters introduces significant challenge due to the high-dimensional combinatorial optimization problem it presents.
We apply our HeO and MCGE to directly train neural networks with ternary value (, 1). Sup- posed that have an input-output dataset D generated by ground-true truth ternary single-layer per-
httpa://biqmac. aau wttp://.aaaaaaaaaaaaaaoaeeeaaeae
http///a....aaaat///aaacccaca.aa www.waatwwwaabbb....... html
https://www. ubc. ca/"hoos/SATLLL benchm. html

---

output dim:
ourpuu dim:
ourput dim: 5
1.0
W
0.8
H40
Mcc
0.6
accuracy
W
0.
Heo {(V. y)...
560
1000
1500 2000
500 1000 1500 2000
500 1000 1500 2000
training set size
training set size
training set size
Figure 5: a, Training networks with ternary-value parameters. b, The weight value accuracy of the HeO (Alg. |1 Appendix) and Monte Carlo gradient estimation (MCGE) with momentum under different sizes of training set (n = 100, m 1,2, 5). We estimate the mean and std from 10 runs.
0.30
0.25
-
q 0.15
q- 0.10
4
4a5
0.9
10
Figure The variable selection of 400-dimensional linear regressions using HeO (Alg. 6 Appendix), Lasso (L,) regression [40] and Lo.5 regression [4T]. We report the accuracy of each algorithm in determining whether each variable should be ignored for prediction and their MSE on the test set.
The mean (dots) and standard deviation (bars) are estimated over 10 runs.
ceptron y = F(v; WGT) = ReLu(WeTv), where ReLu(x) = max æ). WGT E (-1,0, 1,
1) mox is the ground-truth ternary weight parameter, v {-1,0, 13 is the input, and is the model output. We aim to find the ternary configuration W E [-1,0,1]ma
minimizing the loss MSE(W, D) = DD Evyy)E [[T(v; W) - yl1 We generalize our HeO from the binary to the ternary case by representing a ternary variable with two bits, where each element of W can be represented as a function of (see Alg Appendix for details), and the target function defined as
f(s)
Σ (Vy)eD
(T(v;W(()) 2
(16)
D
As shown in Fig. HeO robustly exceeds MCGE under different dataset size |D] and output size m.
Mixed combinatorial optimization. In high-dimensional linear regression, usually only a small fraction of the variables significantly contribute prediction. Identifying and selecting subset of variables with strong predictive power-a process known as variable selection- crucial, as
improves the generalizability and interpretability of the regression model 42 However, direct variable selection is an NP-hard combinatorial optimization mixed with continuous variables [43]
As practical alternative, regularization methods like Lasso algorithm are commonly employed [40]
Supposed a dataset is generated from linear model, which the relation between input v and output y is y 3 + c, where ( is the ground-truth linear coefficient and is independent Gaussian noise with standard deviation a.- Suppose that only small proportion (denoted as (0, 1)) of coordinates in B non-zero. Our goal is to identify these coordinates through an indicator vector 8 E 1- (1 for selection and -1 for non-selection) and estimate these non-zero coefficients. The target function of the problem (1 IR" is the all-one vector)
f(s, ()
1
Σ
(g
8+1
-y
(17)
D
2
(v.))D
8

---

We solve the problem through HeO (Alg. 6. Appendix), where we minimize the loss relative to while slowing varying /3 via its error gradient. After obtaining the indicator s, we conduct ordinary least squares regression on the variables selected by the estimate their coefficients, and treat other variables' coefficients as zero. As shown Fig our HeO outperforms both Lasso regression and the more advanced T...5 regression [4]] in terms of producing more accurate indicators s and achieving lower test prediction errors across various g and o, settings. Importantly, to give the variable selection prediction, our HeO does not need to know the level of q and o, in advance.
covered vertex
graph
examples vertex cover
examples of minimum vertex cover
Figure 7: The illustration of minimum vertex cover.
Table Attributes of graphs and the vertex cover sizes of HeO (Alg
Appendix) and FastVC 44
Graph name
Vertex number
Edge number
FastVC
HcO
tech-RL-caida 190914
607610 1936748
78306
77372 (17)
soc-youtube
495957
149458
148875 (25)
inf-roadNet-PA
1087562
1541514
588306
587401 (104)
inf-roadNet-CA
1957027
2760388
1063352
1061339 (32)
socfb-B-anon
2937612
20959854
338724
312531 (194)
3097165
23667394
421123
387730 (355)
socfb-uci-uni 58790782
92208195
869457
867863 (36)
Constrained binary optimization. Minimum vertex cover (MVC) is class of the constrained combinatorial optimization which has wide applications [45]. as illustrated in Fig. 7 Given an undirected graph with vertex set V and edge &, the MVC to find the minimum subset V. V,
so that for each edge e E &, at least one of its endpoints belongs to Ve. The target function and the constrains are expressed as
8: + 2
8j + subject to (sy(s) = (1 8; (1
((s)
=0, Vi.j.ey EE (18)
2
2
We combine the HcO with penalty function (Alg. Appendix) for solving MVC. We compare HcO with FastVC [44], a powerful MVC heuristic algorithm, on massive real world graph datasets For fair comparison, we keep the time of two algorithms the same for each dataset. As shown in Tab. I. our HeO can find smaller cover sets than that of FastVC.
Discussion
Existing model-based combinatorial optimization approaches encode the solution space via parame- terized distribution with iterative parameter updates [46). In contrast to HeO, which requires only one sample per iteration, they necessitate a large number of samples per iteration. The Gibbs-Wing With- Gradient algorithm 47] uses gradient information for combinatorial optimization but searches in the discrete solution space instead of the continuous one, as did HeO. Denoising diffusion model (DDM) 24 has been applied for solving combinatorial optimization problems
[48). Although the diffusion process in DDMs akin to the heat diffusion in our HeO, DDMs require a substantial data for training and necessitate reversing the diffusion process to generate data that from the target distribu- tion. contrast, HeO needs no training, and is unnecessary to strictly adhering to the monotonic 7 in the optimization process, as under different the function u(T, 0) shares the same optima with that of the original problem h(). This claim is empirically corroborated Fig. S3 Appendix, where HeO applying non-monotonic schedules of 7. still demonstrates superior performance.
http://ettorrrrepstttrycoo

---

Our HeO can be viewed as stochastic Gaussian continuation (GC) method 449 with projection. GC has been applied for non-convex optimization, though it has not yet been used for combinatorial optimization. The optimal convexification of GC [50] underpins potentially theoretical advantages of HeO. One key distinction is that GC typically optimizes each sub-problem (corresponding to u(0, 7x)) at each é up to some criteria, whereas HeO merely performs single-step gradient descent. Also, Eq. (8) corresponds to a variation of the evolution strategy, a robust optimizer for non-diffcrentiable function [5]]. while HcO use a different gradient estimation (Eq. (II), see Appendix for details, Additionally, our HeO is related to randomized smoothing, which has been applied to non-smooth optimization [52] or neural network regularization [53]. The distinctive feature of our HeO is that, across different т, the smoothed function TIT,, 0) retains the optima of the original function h()) (Thm. D.. This distinguishes our HeO from methods based quantum adiabatic theory [9], bifurcation theory [2] and other relaxation strategies [54], in which the optima of the smoothed function can be different from the original one 27 55
This is verified in Fig. S3 Appendix.
The heat equation our HeO can be naturally extended to general parabolic differential equations, given that broad spectrum of them obey the backward uniqueness [56). For example, we can use Ú, 0) = Vel [A V u(T, 0)] to replace the Eq. (6), where A is real positive definite matrix. Prior researches have demonstrated that the optimization procedure can be significantly accelerated by preconditioning [57] or Fisher information matrix (58), implying that choosing a proper matrix A could substantially improve the efficacy of the HeO. Additionally, since Tr does not necessarily have to be monotonic due to the cooperative optimization property of HcO, is feasible to explore various Ti schedules (possibly non-monotonic) to further enhance performance. Moreover, our HcO can be integrated with other existing techniques for combinatorial optimization problems with constraints, such Augmented Lagrangian methods, to achieve better performance (59].
Despite the effectiveness of HeO on the combinatorial optimization problems from different domains we have considered, it has limitations. First, current HeO is inefficient for integer linear programming and routing problems, primarily due to that is cumbersome to encode integer variables through the Bernoulli distribution in our framework. Nevertheless, integrating HeO with other techniques such as advanced Metropolis- Hastoings algorithm [60] may path to broaden the applicability of our methodology to wider range of combinatorial optimization problems. Besides, our HcO allows for further customization by incorporating additional terms that integrate groblem-specifee prior knowledge by hybridizing with other metaheuristic algorithms, allowing for more effective exploration of the configuration space. Second, our HeO can not be theoretically guaranteed for converging to the global minimum. In general, finding the global minimum is not theoretically guaranteed for non-convex optimization problems [61], such as the combinatorial optimization problems studied in this paper. However, it can be demonstrated that the gradient of the target function under heat diffusion satisfies the inequality (22)
C VoulT, 0)1 S VF
where the constant C' depends the dimension. This implies that the target function becomes weakly convex, enabling the finding of global minima and faster convergence under certain conditions [6]
6 Conclusion
In conclusion, grounded in the heat diffusion, we present a framework called heat diffusion optimiza- tion (HeO) to solve various combinatorial optimization problems. The heat diffusion facilitates the propagation of information from distant regions to the solver, expanding its receptive field, which in turn enhances its ability to search for global optima. Demonstrating exceptional performance across various scenarios, our HeO highlights the potential of utilizing heat diffusion to address challenges associated with navigating the solution space of combinatorial optimization.
Acknowledgements
This work was supported by the National Science and Technology Major Project of China (No. 2018AAA0100303) and the Science & Technology Commission of Shanghai Municipality (No. 23JC 23JC1400800).
10

---

References
[1] Francisco Barahona, Martin Grötschel, Michael Jünger, and Gerhard Reinelt. An application of combinatorial optimization to statistical physics and circuit layout design. Operations Research, 36(3):493-513, 1988.
[2] Jun Wang, Tony Jebara, and Shih-Fu Chang. Semi-supervised learning using greedy max-cut. The Journal of Machine Learning Research, 14(1):771-800, 2013.
[3] Chetan Arora, Subhashis Banerjee, Prem Kalra, and SN Maheshwari. An efficient graph cut algorithm for computer vision problems. In Computer Vision-ECCC 2010: 11th European Conference Computer Vision, Heraklion, Crete, Greece, September 5-11, 2010, Proceedings, Purt III 11, pages 552-565. Springer, 2010.
[4] Hayato Sahjjima-Mwesigvaa Christian FA Negre, and Susan M Mniszewski. Graph partitioning using quantum annealing on the d- dave system. In Proceedings of the Second International Workshop on Post Moores Era Supercomputing, pages 22 22-29, 2017.
[5] Florian Neukart, Gabriele Compostella, Christian Seidel, David Von Dollen, Sheir Yarkoni, and Bob Parney. Traffic flow optimization using quantum annealer. Frontiers in ICT, 4:29, 2017.
[6] Román Orús, Samuel Mugel, and Enrique Lizaso. Quantum computing for finance: Overview and prospects. Reviews Physics, 4:100028, 2019.
[7] Rafael Marti and Gerhard Reinelt. Exact and Heuristic Methods in Combinatorial Optimization, volume 175. Springer. 2022.
[8] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, Joshua Lapan, Andrew Lundgren, and Daniel Preda. quantum adiabatic evolution algorithm applied random instances of an np-complete problem. Science, 292(551 41:5516):777-4555 2001.
[9] John A Smolin and Graeme Smith. Classical signature of quantum annealing. Frontiers in physics, 2:52, 2014.
[10] Joseph Bowles, Alexandre Dauphin, Patrick Huembeli, José Martinez, and Antonio Acín. Quadratic unconstrained binary optimization via quantum-inspired annealing. Physical Review Applied, 18(3):034066, 2022.
[11] Mária Erescy-Rave Ravasz and Zoltán Toroczkai. Optimization hardness as transient chaos in an analog approach to constraint satisfaction. Nature Physics, 7(12):966 9700 2011.
[12] Hayato Goto, Kosuke Tatsumura, and Alexander Dixon. Combinatorial optimization by simu- lating adiabatic bifurcations in nonlinear hamiltonian systems. Science advances, 5(4):eaav2372, 2019.
[13] Hayato Goto, Kotaro Endo, Masaru Suzuki, Yoshisato Sakai, Taro Kanao, Yohei Hamakawa, Ryo Hidaka, Masaya Yamasaki, and Kosuke Tatsumura. optimization based on classical mechanics. Science Advances, (16))eabe79533 2021.
combinatorial
[14] Takahiro Inagaki, Yoshitaka Haribara, Koji Igarashi, Tomohiro Sonobe, Shuhei Tamate, Toshi- mori Honjo, Alireza Marandi, Peter L McMahon, Takeshi Umeki. Koji Enbutsu, et al. A coherent ising machine for 2000-node optimization problems. Science, 354(6312):603 2016.
[15] Egor Tiunov, Alexander E Ulanov, and AT Lvovsky. Annealing by simulating the coherent ising machine. Optics express, -222810298, 2019.
[16] Connor Bybee, Denis Kleyko, Dmitri E Nikonov, Amir Khosrowshahi, Bruno A Olshausen. and Friedrich T Sommer. Efficient optimization with higher-order ising machines. Nature Communications, 14(1):6033, 2023.
[17] Martin JA Schuetz, J Kyle Brubaker, and Helmut Katzgraber. Combinatorial optimization with physics-inspiree graph neural networks. Nature Machine Intelligence, 4(4):367-377, 2022.
11

---

[18]
Bernardino Romera-Paredes, Mohammadamin Barekatain, Alexander Novikov, Matej Balog, M Pawan Kumar, Emilien Dupont, Francisco JR Ruiz, Jordan S Ellenberg, Pengming Wang. Omar Fawzi, et al. Mathematical discoveries from program search with large language models. Nature, pages 1-3, 2023.
[19] David Pisinger and Stefan Ropke. Large neighborhood search. Handbook of metaheuristics, pages 99-127, 2019.
[20] Pierre Hansen, Nenad Mladenović, and Jose A Moreno Perez. Variable neighbourhood search: methods and applications. Annals of Operations Research, 175:367-407, 2010.
[21] Haoran Sun, Hanjun Dai, Wei Xia, and Arun Ramamurthy. Path auxiliary proposal for meme in discrete space. In International Conference on Learning Representations, 2021.
[22] Lawrence C Evans. Partial differential equations, volume 19. American Mathematical Society,
2022.
[23] Jean-Michel Ghidaglia. Some backward uniqueness results. Nonlinear Analysis: Theory, Methods & Applications, 10(8):777-790, 1986.
[24] Ling Yang. Zhilong Zhang, Yang Song. Shenda Hong, Runsheng Xu, Yue Zhao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang. Diffusion models: A comprehensive survey of methods and applications. ACM Computing Surveys, 56(4):1 2023.
[25] Endre Boros and Peter Hammer. Pseudo-boolcan optimization. Discrete applied mathematics, 123(1-3):155225) 2002.
[26] Yves Crama and Peter Hammer. Boolean functions: Theory, algorithmss, and applications. Cambridge University Press, 2011.
[27] Bernhard H Korte, Jens Vygen, B Korte, and J Vygen. Combinatorial optimization, volume 1. Springer, 2011.
[28] Shakir Mohamed, Mihaela Rosca, Michael Figurnov, and Andriy Mnih. Monte carlo gradient estimation in machine learning. The Journal of Machine Learning Research, 2181):5183-5244 2020.
[29] David Vernon Widder. The heat equation, volume 67. Academic Press, 1976.
[30] Avradip Mandal, Arnab Roy, Sarvagya Upadhyay, and Hayato Maaijjmma--Meiiwaa.. Com- pressed quadratization of higher order binary optimization problems. In Proceedings of the 17th ACM International Conference on Computing Frontiers, pages 126-131, 2020.
[31] Martin Anthony, Endre Boros, Yves Crama, and Aritanan Gruber. Quadratic reformulations of nonlinear binary optimization problems. Mathematical Programming, 162:115-1444 2017.
[32] George Cybenko. Approximation by superpositions of sigmoidal function. Mathematics of control, signals and systems, 2(4):303-314, 1989.
[33] Michel Gendreau, Yves Potvin, et al. Handbook of metaheuristics, volume 2. Springer, 2010.
[34] John Hopfield and David W Tank. "neural" computation of decisions in optimization problems. Biological cybernetics, 12(3):141-152 1985.
[35] Zhe Wang, Alireza Marandi, Kai Wen, Robert I. Byer, and Yoshihisa Yamamoto. Coherent ising machine based on degenerate optical parametric oscillators. Physical Review A, 88(6))063853, 2013.
[36] Angelika Wiegele. Biq mac library- collection of max-cut and quadratic 0-1 programming instances of medium size. Preprint, 51, 2007.
[37] Holger H Hoos and Thomas Stützle. Satlib: An online resource for research on sat. Sat, 2000:283-292, 2000.
12

---

[38] Andrew Lucas. Ising formulations of many np problems. Frontiers in physies, 2:5, 2014.
[39] Amir Gholami, Schoon Kim, Zhen Dong, Zhewei Yao, Michael W Mahoney, and Kurt Keutzer. survey of quantization methods for efficient neural network inference. In Low-Power Com- puter Vision, pages 291-326. Chapman and Hall/CRC, 2022.
[40] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):267-288, 1996.
[41] Zongben Xu, Hai Zhang, Yao Wang, Xiang Yu Chang. and Yong Liang. L 1/2 regularization. 2010.
Science China Information Sciences,
[42] Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome H Friedman. The elements of statistical learning: data mining. inference, and prediction, volume 2. Springer, 2009.
[43] Xiaoping Li, Yadi Wang, and Rubén Ruiz. A survey on sparse learning models for feature selection. IEEE transactions on cybernetics, 5243):1642-1666, 2020.
[44] Shaowei Cai, Jinkun Lin, and Chuan Luo. Finding a small vertex cover in massive sparse graphs: Construct, local search, and preprocess. Journal of Artificial Intelligence Research, 19:463-494, 2017.
[45] Alexander C Reis, Sean M Halper, Grace E Vezeau, Daniel P Cetnar, Ayaan Hossain, Phillip R Clauer, and Howard M Salis. Simultaneous repression of multiple bacterial genes using nonrepetitive extra-long sgmaa arrays. Nature biotechnology. 37(1 (1::::229411100 2019.
[46] Enlu Zhou and Jiaqiao Hu. Gradient- based adaptive stochastic search for non-differentiable optimization. IEEE Transactions on Automatic Control, 59(7):1818-1888 1832, 2014.
[47] Will Grathwohl, Kevin Swersky, Milad Hashemi, David Duvenaud, and Chris Maddison. Oops i took gradient: Scalable sampling for discrete distributions. In International Conference on Machine Learning, pages 1341-3841. PMLR, 2021.
[48] Zhiqing Sun and Yiming Yang. Difusco: Graph-based diffusion solvers for combinatorial optimization. Advances in neural information processing systems, 2023.
[49] Hossein Mobahi and John Fisher III. A theoretical analysis of optimization by gaussian continuation. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 29, 2015.
[50] Hossein Mobahi and John W Fisher. On the link between gaussian homotopy continuation and convex envelopes. In Energy Minimization Methods Computer Vision and Pattern Recognition: 10th International Conference, EMMCVPR 2015, Hong Kong, China, January 13-16, 2015. Proceedings 10, pages 43-56. Springer, 2015.
[51] Daan Wierstra, Tom Schaul, Tobias Glasmachers, Yi Sun, Jan Peters, and Jürgen Schmidhuber. Natural evolution strategies. Journal of Machine Learning Research, 15(27):949-99) 9800, 2014.
[52] John C Duchi, Peter I. Bartlett, and Martin J Wainwright. Randomized smoothing for stochastic optimization. SIAM Journal on Optimization, 22(2):674-701, 2012.
[53] Mo Zhou, Tianyi Liu, Yan Li, Dachao Lin, Enlu Zhou, and Tuo Zhao. Toward understanding the importance of noise in training neural networks. In International Conference Machine Learning, pages 7594-7602. PMLR, 2019.
[54] Nikolaos Karalias and Andreas Loukas. Erdos goes neural: an unsupervised learning framework for combinatorial optimization on graphs. Advances in Neural Information Processing Systems, 33:6659-6672, 2020.
[55] Haoyu Peter Wang, Nan Wu, Hang Yang, Cong Hao, and Pan Li. Unsupervised learning for combinatorial optimization with principled objective relaxation. Advances in Neural Information Processing Systems, 35:31444-31458, 2022.
[56) Jie Wu and Liqun Zhang. Backward uniqueness for general parabolic operators the whole space. Calculus of Variations and Partial Differential Equations, 18:1-19, 2019.
13

---

[57] Hengyuan Ma, Li Zhang, Xiatian Zhu, and Jianfeng Feng. Accelerating score-based generative models with preconditioned diffusion sampling. European Conference on Computer Vision, pages 1-16. Springer, 2022.
[58] Daan Wierstra, Tom Schaul, Tobias Glasmachers, Yi Sun, Jan Peters, and Jürgen Schmidhuber. Natural evolution strategies. The Journal of Machine Learning Research, 15(1):949 -800, 2014.
[59] Ernesto G Birgin and José Mario Martínez. Practical augmented Lagrangian methods for constrained optimization. SIAM, 2014.
[60] Haoran Sun, Katayoon Goshvadi, Azade Nova, Dale Schuurmans, and Hanjun Dai. Revisiting sampling for combinatorial optimization. In International Conference on Machine Learning, pages 32859-32874. PMLR, 2023.
[61] Feng-Yr Liao, Lijun Ding, and Yang Zheng. Error bounds, pl condition, and quadratic growth for weakly convex functions, and linear convergences of proximal point methods. In 6th Annual Learning for Dynamics Control Conference, pages 993-1005. PMLR, 2024.
[62] Felipe Atenas, Claudia Sagastizábal, Paulo JS Silva, and Mikhail Solodov. A unified analysis of descent sequences in weakly convex optimization, including convergence rates for bundle methods. SIAM Journal on Optimization, 33(1):89-115, 2023.
[63] Paul Bertens and Seong-Whan Lee. Network of evolvable neural units: Evolving to learn at synaptic level. arXiv preprint, 2019.
[64] Juntao Wang. Daniel Ebler, KY Michael Wong, David Shui Wing Hui, and Jie Sun. Bifurcation behaviors shape how continuous physical dynamics solves discrete ising optimization. Nature Communications, 14(1):2510, 2023.
14

---

Appendix
A Proof of theorems
A.1 Proof of Thm.
To prove the Thm. we recall the backward uniqueness of the heat equation 23 which asserts that the initial state of a heat equation can be uniquely determined by its state a time point T under mild conditions, as shown in the following theorem.
Theorem S4. Given two bounded function h1 (x), b2 (x) with domain on R Denote 41 (T, x) and (T x) the solutions the heat equation (Eq. (6)) with initial condition un (0, x) h1 (x) and
(2 x) hx(x) respectively. If there exists a r > 0, such that
(T, ()) E
(S1)
we have h1 = h22-
Proof of the Thm.
Proof. We first show that the global minima of u(T, 0) ] is also the global minima of h()). The corner- stone of the proof is Thm. [S4] To utilize the backward uniqueness, we consider a reparameterizatioo of p(s}0) by introducing a random variable x ~ Unif(O, 1)". It is easy to see that sgn(O; x..) obeys to be 1. Replacing 8; by x, in Eq. (2)
Bernoulli distribution with probability 0, to be and we have new expression of (0)))
h(0)
x))], x (0,
(S2)
Since the solution to the heat equation can be represented by the heat kernel, we have 22
" 0) (h(0 + 22T2)], p(z) = I).
(S3)
Therefore, we have
u(T, 0) p(z (x)) (sgn /2Tz - p(xx) (ppp)) I [sgn(( (x + V2T2))]]-
(S4)
Denote u(T, x: 0) = Ex(x) Of (sgu(0 - +x+ v2Tz)))[, x E (0, 1)". Noticed that u(T, x; 0) is the solution of the following unbounded heat equation respect to the time and location restricted on the region x E 10,
{
( x; = Axu(T, x: 0), u(TT,x
T T>0, T T= 0,
x x R" x E R
(S5)
f(sgn(0 - x)),
hence the latter can be considered as extension of the former. Since u(T, 0) is analytic respect to x E 1)" for > 0, this extension is unique. Therefore, the value of u(T, x; 0) on (0, 1)", > uniquely determines the solution of Eq. (S5) Denote as
+00, sf +1
(S6)
0o, sf 1.
Then u(T,x:0) =-, for any and R", and we have
u(T, 0)
f
(S7)
Noticed that for T we have
u(r, 0) Eppx [M(T, x; 0)] Ex [B (Sgu((O (x V2T2)))]] < (Epc)) [s]]] [[]]=]
(S8)
and the equality is true if and only u(T, x; 0) = f is true for x E R", Therefore, if θ is the one of the minimas of u(T, 0), and we have u(T, 0) > f*. Similarly, since
u(r, 0) Ep(x) (u(T, x: 0)] (Bp (sga( (x v2rz)))Il < Ep(x) (5p())11) f
(S9)

---

we also have u(T,0) S f", hence
= u(T, 0*), = x R".
(S10)
Due the backward uniqueness of the heat equation, we have
u(0, x: 0) x; 0*), x KER,
(S11)
that is
h(0) h(0)) f
(S12)
As result, f is the of minimas of h()). Conversely, using Eq. i is one of minimas of h()), it is also one of minimas of u(T, 0).
(S7 it obviously to that if
A.2
Proof of Thm.
Recall Eq.
S3 and use the definition of h())) (see Sec. 2 of the main paper), we have
u(T, 0) Ep(x) [ p(2)0++222)) [s])]]
(S13)
To construct low-variance estimation, instead of directly using the Monte Carlo gradient estimation by sampling from p(z) and p(s|0 + 222z) for gradient estimation, we manage to integrate out the stochasticity respect z. Use the reparameterization Eq. (S2) and Eq. (S3). we have
u(T 0) Ep(x) (pt O v2zz x)]]]
(S14)
Now we can calculate the inner term Ep (sgn(( + /2Tz - x))]. Due to the assumption, the target function (() can written as a K-order multilinear polynomial of
((8) = ao + Σ
an 81g +
Σ
81, +
Σ
81 81,
81
(S15)
(Kii....g
Six
i <in
Integrating respect to each dimension of the Gaussian integral, we have [50
(fte))((sgn(( + /2Tz x))]
ag + Σ a1.i Si
Σ
a2,i 5 S1
Σ
03,i
S1,
(S16)
aKi Bi 8ix
where
0: Tj
3, = px(( (sgn(0; + 2 x.)] erf(
(S17)
/
where erf(-) is the error function. Therefore, have
u(T, 0) [f x VT })}
(S18)
where e-) is the element-wise error function.
A.3 Proof of Thm.
Proof. Define the square loss as
c()) = (h(0) -
(S19)
According to the definition of h()),, we have
e(0) (f - -x)) (0 -x)]]² < Ep(x) (f(ggn( x)) f x))
(S20)

---

Define the error function
r(T, 0) = u(T, 0) u(T, x; 0*).
(S21)
Then the enror function satisfies the following heat equation
{
(i,r(T,x:0) r(0, 0)
Vxt(T, 0)
(S22)
fssgu(O 0 - -x) f(sgn(O* x))
Define the energy function of the error function r(T, x; 0) as
/ xT,xx;0)(xx))
(S23)
JR"
Then applying the heat equation and the integration by parts, we have
E(T:0) 2 / x)]] p(x)dx. dr
(S24)
Hence we have for 0 < T2
+)+ 7. 1/( ||vr(T, 0)11 p(x)dxxT.
B( 0) E( 0) +
(S25)
Use the Harnack's inequality [22] we have
||r[(, 0) S 0(),,,
+ n (T; 0),
(S26)
2r
combine with Eq. S25
we have
E(T; dr.
E(T1;)) S E(T2; 0) + 2
(S27)
Using the Minkowski inequality on the measure p(x), we have
((0) - h(00) = e2)2((0) < ( { (f(sgn(0 x)) - x(T;;x 1/2
(u(TT:x;0) u(1:: 0)) p(x))xx 1/2
+
R"
(S(sgn(O" x)) - ( 1/2
+
(S28)
R"
= (f (f(sgu(0 x)) - u(TT: 0)) (x)x)xx 1/2
f (f(sgn(0" x)) - 0")xx 1/2 (1:)
R"
Recall the continuity of the heat equation:
lim RR
(u(T, 0) fsgn -x)))=(()xxxxxxxxx))) 0.
(S29)
Therefore, given 0, there exists a Ti 0, such that
(ff (J(sgn(0- x)) ( 1/2
(S30)
+ ( / (f(sgn(O* - x)) - (z) 1/2 <E.
Recall Eq. S27). we then have the error control for e(0):
e/22 < (21(((1:) +es (E(T2: 0) ( E(T; dr) 1/2
+e.
(S31)
Noticed that
E(T: < ($ - E') u(T,x; (")((T (0,0))--66666 O((,)")),
(S32)
where j = maxx f(s). and we prove the theorem.
3

---

40
time per iteration (ms)
20
15
10
01
0.0
0.5
1.0
1.5
2.0
2.5
3.0
dimension
109
Figure S1: The time cost per iteration (ms) of the HeO framework increases linearly with the dimensionality of the problem. We present the results averaged over five tests, with error bars representing three standard deviations.
gka
rudy
40
Figure S2: Mean result target function values (lower is better) our HeO on three datasets (Ising, gka, and rudy) under various step sizes (7) and iteration counts (7). The star denotes the settings used Figure of the main paper.
B Complexity analysis
We empirically estimate the time cost per iteration of our HeO for problems of different dimension Fig. SI
C Parameter sensitivity analysis
We report the performance of our HeO wide range of parameter settings in Fig. $2

---

D Relation to evolution strategy
We show that the update in Eq. (8) is equivalent variation of the evolution strategy (ES) which is a robust and powerful algorithm for solving box optimization problem 511 In fact, we have
Vou(T, 0) Vo ppx)
[h(0 + v2Tz)]
Vo
exp(
1|z - (()d)
V ATT
47
exp(
1z 01 013))((2)( 00
2T
V4
4T
(S33)
1
exp(-
h2 + 22T2)zdz
2T
/2T"
px)) [h(0 2Tz)zz]
The random vector z corresponds to the stochastic mutation and h(O v2rz) corresponds to the fitness Alg. in [63]. In ES, the standard deviation 2T is fixed in general, while our HeO v 2T. is varying across time. As shown in Thm. 3 the varying T: in our HeO offers theoretical benefit that controls the upper bound this benefit. Another difference between our HeO and ES is that we integral out and using Eq. to estimate the gradient, while ES the gradient is estimated based on sampling z.
optimization result. In contrast, constant t in ES does not provide
E Relation to denoising diffusion models
Our approach, while bearing similarities to the DDM highly regarded and extensively utilized artificial generative model [24] that relies on the reverse diffusion process for data generation differs in key aspects. The DDM necessitates reversing the diffusion process to generate data that from the target distribution. sequence in the optimization process, as under different the function u(T. (T. 0 shares the same optima with that of the original problem h(0). This claim corroborated in Fig. S3 as shown below, where HeO applying non-monotonic schedules of still demonstrates superior performance. Hence, is possible to explore diverse T, schedules to further performance enhancement.
contrast, is unnecessary for our HeO strictly adhering to the reverse time
F Cooperative optimization
30000
25000
LGA
20000
bSB
d58
15000-
aSB
10000-
SIM
5000
Best known
Figure S3: Verifying the cooperative optimization mechanism of HeO. The best cut value over 10 runs for each algorithm on the K-2000 problem [14] when the control parameters are randomly perturbed by different random perturbation level 5. The red dash line is the best cut value ever find.
Our study suggests that HeO exhibits a distinet cooperative optimization mechanism, setting it apart from current methodologies. Specifically, HeO benefits from the fact that target functions u(T, O) share the same optima as the original problem h(0) for any > 0. This characteristic allows the solver to transition between different 7 values during the optimization process, eliminating the necessity for monotonic schedule. In contrast, traditional methods such as those based on quantum adiabatic theory or bifurcation theory require a linear increase of control parameter a from 0 to 1. This parameter is analogous to 7; in the HeO framework.
5

---

To empirically verify the above claim, we introduce a random perturbation to the schedule in Alg. 1, rendering it non-moontooni:: 7, T. where e uniformly distributed on [1 6, 1++ 6] with controlling the amplitude of the perturbation. For other methods based on quantum adiabatic theory or bifurcation theory, we correspondingly introduce the perturbation as 0₂ = clamp(c,a, 0, 1). If an algorithm contains cooperative optimization mechanism, it still works well even when the control parameter is not monotonic, as optimizing the transformed problems under different control parameters cooperatively contributes optimizing the original problem. As shown in Fig. [3] the performance of other methods are all dramatically deteriorated. In contrast, HeO shows no substantial decline in performance, corroborating that HeO employs cooperative optimization mechanism.
G Implementation details
All the experiments are conducted single NVIDIA RTX-3090 GPU (24GB) and an Intel(R) Xeon(R) Gold 6226R CPU (2.90GHz). For convenience, we denote VT, = 0z.
Monte Carlo gradient estimation for combinatorial optimization. Based on Eq. (5) we construct a combinatorial optimization algorithm using the Monte Carlo gradient estimation (MCGE), as shown in Alg. 2 Noticed that we clamp the parameter 0 in the [0. 1] for numerical stability; additionally, we binarize the rr to obtain the optimized binary configuration ST the end.
Algorithm 2 Monte Carlo gradient estimation for combinatorial optimization (MCGE)
Input: target function f-), step size y, sample number M, iteration number T initialize elements of Og as 0.5 for t = 0 to T - 1 do
sample (m)) m..idd p(s 0c); m m= I.....,,M M
g: < M Σm (m(m)) Vo, logp(s(m)(0))
0,-1 t clamp O2 78., 0,1)
end for
f- sgn(Or 0r - 0.5) Output: binary configuration ST
Gradient descent with momentum. We provide HeO with momentum in Alg. 3
Algorithm 3 Heat diffusion optimization (HeO) with momentum
Input: target function f(-), step size 7, momentum k, schedule {ar) iteration number T' set g-1
initialize elements of 0, as 0.5 for t T 1 do sample x, ~ Unif (0, 1]"d
w. Vo, f (erf( (==xe)
g: < kg: 1 1 We
0,11 - Projz 0 + gx)
end for
ST sgu(Or 0.5) Output: binary configuration ST
Toy example. We set the momentum x = 0.9999, learning rate 7 = and iteration number T = 5000 for HeO (Alg. 3. For MCGE without momentum, we set T 50000, set =le-6, momentum k = 0, and M 10. For MCGE with momentum, we set momentum as 0.9.
Max-cut problem. For solving the max-cut problems from the Biq Mac Library [6], we set the steps T = 5000 for all the algorithms. For HeO, we set y and o, linearly decreases from 1 to for HeO, and we set momentum as zero. For LQA and SIM-CIM, we use the setting in [10]. For bSB, dSB, aSB, and CIM, we apply the settings algorithm alg and instant the relative loss is calculated as |C the lowest output of the algorithm alg on the instance i over tries, and Ciii the lowest output of all 7 the algorithm the instance i. For each test, we estimate the mean and std from 10 runs.
664
To reduce the fluctuations of the results, for each
where Chaag is

---

3-SAT problem. For Boolean 3-satisfiability (3-SAT) problem, we set the momentum x = 0.9999, T' = 5000, 7 2, and linearly decreases from v2 to 0 for HeO, According to the empirical finding that high-order loss function leads to better results []6] we include higher order terms in the target function. However, since
1 4;
(S34)
2
for any c; E (-1,1}, directly introducing higher order term in the f(s) is useless. Instead, we adjust the gradient from
0: X1 1))
H
3
ch. (erf( 0: -x)x)))))
Verf(erf(
Ve,
)) )
(S35)
2
H 3
0: 0x=x1)) )x) ))*} Vo. Σ1 2 Ch. (erf( A=1==
(S36)
dt
We consider the 3-SAT problems with various number of variables from the SATLIB 37 For each number of variables in the dataset, we consider the first 100 instance. apply the same configuration of that in [16] for both 3-order solver and 2-order oscillation Ising machine solver. The energy gap of the 2-order solver is set as 1. For each test, we estimate the mean and std from 100 runs.
Algorithm 4 Heat diffusion optimization (HeO) for 3-SAT problem
Input: adjusted target function f(s) = Σm= II:
1-ca,#h,
(Eq. S36 ), step size y, momen-
2
tum K, a schedule { iteration number T
initialize elements of 0 as 0.5, set g-1 0 for t 7 do
sample x, Unif 1Jmd
Va, / (erl( [ (=xx))
g: < kg + W de
0.-1 t Projz 0, 7gr)
end for
ST sgn(Or- 0.5) Output: binary configuration ST
Ternary-value neural network learning. We represent a temary variable E (,1,0,1) as s (sb, + sb.2) with two bits $6,2 E {-1,1). In this way, each element of W can be represented as function of s [mxxxx2 We denote this relation matrix-value function W = (s))
(s) = (stj. + s15,2), i= i, m, .J 1,..,,
(S37)
Based on the above encoding procedure, we design the training algorithm for based on HeO in Alg 5 The input of the dataset D is generated from the uniform distribution (-1, 1)". For HeO, set T 10000, 7 0.5, 6 0.999, and o, linearly decreasing from /2 to 0. For MCGE, set T = 10000, 7 = le 7, M = 10, and x = 0.9999. We empirically find that MCGE need high sampling number (M and low learning rate (ry) for stability, while this is not the case for HeO. For each test, we estimate the mean and std from 10 runs.
Variable selection problem. We construct algorithm for variable selection problem based on HeO as shown in Alg. 6, where the function f(s, (3) is defined in Eq.
We randomly generate di0-dimensional datasets with 1000 training samples. The input v is sampled from a standard Gaussian distribution. The element the ground-truth coefficient 3 uniformly distributed on 2-2, -1] U [1, 2]. and each element haves 1 1-1 probability of being set zero and thus should be ignored for the prediction. We apply a five-fold cross-validation for all of methods. For our HcO, we set T 2000 and 1, 0.999. We generate ensemble of indicators of size 100.
7

---

Algorithm 5 Heat diffusion optimization (HeO) for training ternary-value neural network
Input: dataset D, Step size 7, momentum K, schedule {oe), iteration number T
initialize elements of 0g as 1, initialize elements of 30 as O, set g = 0 for /. = 0 to 7 do
sample X Unif (0, 1)
Wr - (erf(
(Eq. S377
MSE D Eivy)ED T(v; W y12
W: - V,, MSE
g: f kg:-1 -1 W
0c-1 f Projz (, g:)
end for
ST = sgn (0r - 0.5)
Output: Wr ((r)
Algorithm 6 Heat diffusion optimization (HeO) for linear regression variable selection
Input: dataset D, step size momentum K, schedule {or), iteration number T
initialize elements of 0g as 1, initialize elements of 30 as O, set g = =g" 0.
for t = 0 to 7 1 do
sample x 00 11Td
w₂" V3, / [(rff(
1,
g: k Tw?
9
G - B, gi
w + Vo, Serf( 0e-Xx ).
gi kg! + yw,
0-11 Projz O0 - 7g()
end for
ST f sgn(0; 0.5)
Output: ST
For each s in the ensemble, we fit a linear model by implementing an OLS on the non-zero variables indicated by and calculate the average MSE loss of the linear model on the cross-validation sets. We then select the linear model with lowest MSE on the validate sets as the output linear model. For Lasso and To.55 regression, we follow the implementation in [4]] with 10 iterations. the regularization parameter is selected by cross-validation from (0.05, 0.1, 0.2, 0.5, 1, 2, 5). For each test, we estimate the mean and std from 10 runs.
Table S2: The attributes of the real world graphs and the parameter settings of HeO.
graph name
VV
F|
T
tech-RL-caida
190914
607610
200
2.5
linearly from to 2.5
soc-youte soc-youtube
495957
1936748
200
2.5
lincarly from to 2.5
inf-roadNet-PA,
1087562
1541514
200
2.5
lincarly from 0 to 7.5
inf-roadNet-CA
1957027
2760388
200
linearly from 0 to 7.5
linearly from v2 0
boff-B-anoo
2937612
20959854
50
2.5
linearly from to 5
soefb-A-anon
3097165
23667394
50
2.5
linearly from to
soefb-uci-uni
58790782
92208195
50
2.5
linearly from to
Minimum vertex cover problem.
For constrained binary optimization
min f(s), sE(1,1
(S38)
(s) < 0, k = 1....,A K,
(S39)
8

---

we put the constrains as the penalty function with coefficient À into the target function
K
(s(s) f(s)
ghe (s),
(S40)
and the corresponding algorithm is shown in Alg.
Algorithm 7 Heat diffusion optimization (HeO) for constrained binary optimization
target function with penalty Jx, step size 7, a schedule (ax), penalty coefficients schedule (A+). iteration number T initialize elements of 0g as 0.5, set get 0. for t = 0 to T 1 do sample x from Onif(o, 1]"
Input:
W; f Ve. Sx, (x=X))
g: < kg. WW
0,-1 - Projz 0₁ 7gv)
end for
ST f- sgn(0s (r 0.5)
Output: ST
Algorithm 8 Refinement of the result of MVC
Input: the result of HeO ST
for i = 1 to n. do
set ST, as 0 if Sy is still a vertex cover
end for Output: ST
We implement the HeO single NVIDIA RTX 3090 GPU for all the minimum vertex cover (MVC) experiments. Let s be the configuration to be optimized, in which s; is I if we select i-th vertex into Vo, otherwise we do select i-vertex into V.. The target function be minimize the size of Ve: f(s) = >i i=
1 and the constrains are
8, 1
8; + (1
Pig (s) = (1
=0
Vi,jey E
(S41)
where ey represent the edge connecting the i and j-th vertices. We construct the target function (s) - (() + > Enyer (s). The term with the positive factor penalizes vector when there are uncovered edges. After the HeO outputs the result ST. we empirically find that its subset may also form a vertex cover for the graph G, so we implement the following refinement on the result sy, as shown in Alg. 8 We report the vertex number, edge number and settings of HcO in Tab. S2 For VastVC, we follow the settings in [4] and use its codebase, and set the cut-off time as the same as the time cost of HeO. For each test, we estimate the mean and std from 10 runs.

---

NeurIPS Paper Checklist
1. Claims
Question:
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We propose a framework for solving general combinatorial optimization problems based heat diffusion called HeO. We provide theoretical explanation and exhaustive experiments to demonstrate the superior performance our framework on combinatorial optimization problems across different fields.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation as long as is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We discuss limitations in the last paragraph of the Sec. 5.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "I "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
10

---

Answer: [Yes]
Justification:
We provide detailed proofs in the Appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: We provide information of reproduction in Implementation details section of Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e-g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
The codes are provided the supplementary files.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
//nips //iiss.cc/ppll wttpcccpppppppppeeeeeeeeeeeeeeee httaa//gggiiee//eeeuebiisisssss /Cc////oddiiisssseeeeeee ssiiiiiiccyyyyyyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification:
We provide experiment details in Implementation details section of Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We provide error bars for all experiments in Fig. 3-6 and reports the std in Tab.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for
example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
12

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative rates).
show in tables or
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We provide experiment details in Implementation details section of Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't It make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer:
[Yes]
Justification:
The research conducted in the paper conform with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- cration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: The research on combinatorial optimization has no immediate social impact or potential harms.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
13

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper does not contain pretrained language models, image generators, scraped datasets.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
The URLs of the datasets are provided in the main paper.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
14

---

If assets are released, the license, copyright information, and terms of use in the thcode. has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, batasets.papersstt paperswi
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors encouraged to reach out the asset's s creators.
13. New Assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: The released code is well-documented, and the documentation can be found in the README.md file within the codebase.
Guidelines:
The answer NA means that the paper does not release assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
15

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16


Reviews:
[{"review_id": "ecT4DMTnwa", "reviewer": "Reviewer_NRJw", "summary": "The paper presents a novel framework for solving combinatorial optimization problems using a concept termed \"Heat diffusion optimization (HeO).\" The approach diverges from traditional methods by utilizing heat diffusion to enhance information propagation within the solution space, allowing for more efficient problem-solving.", "questions": "1.\tOn Page 4, Line 129, the method for determining the value of K in the multilinear polynomial is not clear. Could you elaborate on how K is selected in different scenarios?\n2.\tFigure 2 shows that the energy does not monotonically decrease during optimization. Why the energy in Fig. 2 is not monotonically decreased?\n3.\tWhat are the algorithmic complexity and computational time of the proposed heat diffusion optimization method?\n4.\tIs it possible to apply your heat diffusion framework to other types of combinatorial optimization problems, such as those found in operational research?", "limitations": "1.\tOn Page 3, line 105, the definitions of \\Delta_\\theta and the function u are not given.\n2.\tThe ability of global search is not theoretically analyzed.", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1.\tThe introduction of heat diffusion as a mechanism to aid in combinatorial optimization is novel and thoughtfully developed.\n2.\tComprehensive experiments across different optimization problems illustrate the method's effectiveness and superiority over traditional approaches.\n3.\tThe methodology is tested on a wide range of problems, showing its versatility and potential for broader application in real-world scenarios.", "weaknesses": "1.\tThe paper lacks a detailed discussion on the scalability of the method, especially in very large-dimensional spaces, which are common in real-world applications.\n2.\tMore baselines are needed.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 1. Could you elaborate on how K (Page 4, Line 129) is selected in different scenarios?\n\n\nThank you for your valuable feedback. We will clarify this point after revision. The value of K is directly determined by the target function $f(\\mathbf{s})$ to be optimized. For example:\n\nFor minimum vertex cover problem, where $f(\\mathbf{s})$ is a linear function (Eq. 18), we have K=1.\n\nFor Quadratic unconstrained binary optimization, where $f(\\mathbf{s})$ is a quadratic function (Eq. 14), we have K=2.\n\nFor 3-satisfiability problem, where $f(\\mathbf{s})$ is a cubic function (Eq. 15), we have K=3.\n\n> 2. Why does the energy in Fig. 2 not monotonically decrease?\n\nIn simulated annealing, non-monotonicity arises due to the stochastic nature of the annealing process, which occasionally permits transitions from lower to higher energy states. \n\nIn MCGE and our proposed HeO method, the non-monotonic behavior is a result of the stochasticity of the gradient estimate (Eq. 5 and Line 4-5 of Alg. 1 in the main paper).\n\n> 3. The algorithmic complexity, computational time, and the scalability of the HeO, especially in very large-dimensional spaces.\n\nThank you for pointing this out. We will include a discussion about the computational complexity in the main paper. The complexity of our algorithm is relatively low. At each step, the most computationally intensive operation is calculating the gradient (Line 5 of Alg. 1 in the main paper). This can sometimes be expressed explicitly or be efficiently computed using automatic differentiation tools like PyTorch's autograd. The overall computational time is primarily dependent on the number of iterations $T$. As demonstrated in Fig. 1 of the ***PDF file of the global author rebuttal***, the time cost per iteration of our methods increases linearly with the problem dimension, with a small constant coefficient. This confirms the efficiency of our method. We will include this result in the Supplementary Information after revision.\n\n> 4. Is it possible to apply your heat diffusion framework to other types of combinatorial optimization problems?\n\nYes, our HeO framework is versatile and can be applied to various types of combinatorial optimization problems. Since QUBO, 3-SAT, and Minimum Vertex Cover are NP-complete problems, any NP-hard combinatorial optimization problem can theoretically be encoded into these forms and then solved using HeO. For instance, our framework can be used to solve problems like the multidimensional knapsack problem and graph coloring.\n\n> 5. More baselines are needed.\n\nWe have included representative baselines across various combinatorial optimization problems, including satisfiability (e.g., 3-SAT), graph theory (e.g., max-cut, minimum vertex cover), neural network training (e.g., ternary network training), and statistics (e.g., variable selection for linear regression). In our comparisons, the HeO framework performed better against advanced methods, including state-of-the-art higher-order Ising machines proposed last year [2], as well as coherent Ising machines (CIMs) and its variants, which are widely used in industry [3]. Our goal is not to claim that HeO surpasses specialized methods for specific combinatorial optimization tasks, as metaheuristics can be tailored to achieve state-of-the-art results for particular problems [4]. Instead, we aim to highlight the generality of HeO: *it offers a different theoretical approach while performing competitively across a wide range of combinatorial optimization problems without requiring specialized design* for each case. This makes HeO a promising framework for addressing a broad spectrum of combinatorial optimization challenges.\n\n> 6. On Page 3, line 105, the definitions of $\\Delta_\\theta$ and the function $u$ are not given.\n\nThank you for pointing this out. We will clarify this after revision. The function $u$ is the solution to the heat equation as defined in Eq. (6), and $\\Delta$ is the Laplace operator $\\Delta_{\\mathbf{x}} f(\\mathbf{x}) = \\sum_{i=1}^{n} \\frac{\\partial^2 f}{\\partial x_i^2}$.\n\n> 7. The ability of global search is not theoretically analyzed.\n\nThank you for pointing this out. We will include a theoretical discussion after revision. In general, finding the global minimum is not theoretically guaranteed for non-convex optimization problems [5], such as the combinatorial optimization problems studied in this paper. However, it can be demonstrated that the gradient of the target function under heat diffusion satisfies the inequality [6]: \n\n\\begin{align}\n\\left\\| \\bigtriangledown_{{\\theta}}u(\\tau,{\\theta} )\\right\\| \\leq \\frac{C}{\\sqrt{\\tau}},\n\\end{align} \n\nwhere the constant $C$ depends on the dimension. This implies that the target function becomes weakly convex, enabling the finding of global minima and faster convergence under certain conditions [7].\n\n[1] Korte B H, Vygen J, Korte B, et al. Combinatorial optimization[M]. Berlin: Springer, 2011.\n\n[2] Bybee C, Kleyko D, Nikonov D E, et al. Efficient optimization with higher-order Ising machines[J]. Nature Communications, 2023, 14(1): 6033.\n\n[3] Wang J, Ebler D, Wong K Y M, et al. Bifurcation behaviors shape how continuous physical dynamics solves discrete Ising optimization[J]. Nature Communications, 2023, 14(1): 2510.\n\n[4] Glover, Fred W., and Gary A. Kochenberger, eds. Handbook of metaheuristics. Vol. 57. Springer Science & Business Media, 2003.\n\n[5] Liao F Y, Ding L, Zheng Y. Error bounds, PL condition, and quadratic growth for weakly convex functions, and linear convergences of proximal point methods[C]//6th Annual Learning for Dynamics & Control Conference. PMLR, 2024: 993-1005.\n\n[6] Evans L C. Partial differential equations[M]. American Mathematical Society, 2022.\n\n[7] Atenas F, Sagastizábal C, Silva P J S, et al. A unified analysis of descent sequences in weakly convex optimization, including convergence rates for bundle methods[J]. SIAM Journal on Optimization, 2023, 33(1): 89-115.", "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": "1VstfFOoLu", "reviewer": "Reviewer_c1oa", "summary": "This paper aims to improve the efficiency of existing combinatorial optimization methods via heat diffusion. The author have made a thorough analysis over the existing problems and propose the heat diffusion method HOE for general combinatorial optimization problems. The empirical evaluation verifies its advantage over various combinatorial optimization problems.", "questions": "See the weaknesses part", "limitations": "Yes", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The authors make a comprehensive analysis of the problems of existing methods, the proposed method is quite novel, with enough insights to the future research on the combinatorial optimization.\n2. The proposed method is both theoretically supported and empirically justified.\n3. The empirical evaluation spans a various of combinatorial optimization problems, and thorough analyses are presented.\n4. The whole paper is well-written", "weaknesses": "1. In section 2, the analyses are only focused on the methods doing the gradient descent over the relaxed variables. It is unclear to me how the conclusions are generalized to the method like large neighborhood search, variable neighborhood search and path auxiliary sampling (as mentioned in the introduction part).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> In section 2, the analyses are only focused on the methods doing the gradient descent over the relaxed variables. It is unclear to me how the conclusions are generalized to the method like large neighborhood search, variable neighborhood search and path auxiliary sampling.\n\nGood question. We acknowledge that generalizing the analysis from single-step search to multi-step search scenarios is challenging. In this paper, we focus on single-step search because, while multi-step search can improve the chances of escaping local minima and finding better solutions in general, it also increases computational costs, requires more careful design of search rules, and carries the risk of backtracking problems [1]. Integrating our HeO framework with multi-step search methods could be a valuable direction for future research, both theoretically and practically.\n\n[1] Sun H, Dai H, Xia W, et al. Path auxiliary proposal for MCMC in discrete space[C]//International Conference on Learning Representations. 2021.", "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": "yWcjLchS35", "reviewer": "Reviewer_uVQ4", "summary": "This paper proposes the Heat Diffusion Optimization (HeO) method, which leverages thermodynamic principles to enhance combinatorial optimization (CO) problems. Specifically, it integrates heat diffusion equations into gradient-based optimization to improve efficiency and help escape local minima.", "questions": "None.", "limitations": "No concerns.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "1. The paper is generally well-written, with a good balance of examples, explanations, and discussions.\n2. The idea of using heat diffusion to propagate information across the solution space is novel and interesting for solving CO problems.\n3. The proposed method is validated on different types of CO problems with varying scales.", "weaknesses": "1. While the work introduces a new approach for gradient-based combinatorial optimization, it inherits the limitations of gradient-based methods. It may struggle with complex CO problems, such as routing problems, as discussed by the authors.\n\n2. The paper could benefit from a deeper theoretical and empirical analysis of the HeO algorithm. For example, a detailed analysis of the convergence properties and computational complexity of the algorithm is needed. Also, summarizing and providing recommendations on parameter sensitivity and selection would be valuable.\n\n3. The paper demonstrates improved performance over several classic solvers like MDGE, SA, and LQA. However, it does not provide enough evidence or discuss the potential of HeO to achieve state-of-the-art performance on the studied CO problems.\n\n4. The submission seems to be not incomplete as there lack of appendices or supplementary materials that verify the soundness of the theorems presented in Section 3.\n\nOverall, I feel like this paper is proposing a very interesting new method for CO with potential. It would benefit greatly from a major revision that includes more details, proofs, and broader experiments.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 1. The paper could benefit from a deeper theoretical and empirical analysis of the HeO algorithm. For example, a detailed analysis of the convergence properties and computational complexity of the algorithm is needed. Also, summarizing and providing recommendations on parameter sensitivity and selection would be valuable.\n\n\nThank you for the valuable suggestions.\n\nFor **convergence analysis**, we will include a theoretical discussion in the main paper. In general, finding the global minimum is not theoretically guaranteed for non-convex optimization problems [1], such as the combinatorial optimization problems studied in this paper. However, it can be demonstrated that the gradient of the target function under heat diffusion satisfies the inequality [2]: \n\n\\begin{align}\n\\left\\| \\bigtriangledown_{{\\theta}}u(\\tau,{\\theta} )\\right\\| \\leq \\frac{C}{\\sqrt{\\tau}}, \n\\end{align}\n\nwhere $C$ is a constant depended on the dimension and $||f||_{\\infty}$. This suggests that the target function becomes weakly convex, which has been shown to enable the discovery of global minima and achieve faster convergence rates under certain conditions [3].\n\n\nFor **complexity analysis**, we will add a discussion on computational complexity in the main paper. The complexity of our algorithm is relatively low. At each step, the most computationally intensive operation is calculating the gradient of the target function $h(\\theta)$. This can sometimes be expressed explicitly or be efficiently computed using automatic differentiation tools like PyTorch's autograd. The overall computational time is primarily dependent on the number of iterations $T$. As demonstrated in Figure 1 in the **PDF file of the global author rebuttal**, the time cost per iteration of our methods increases linearly with the problem dimension, with a small constant coefficient. This confirms the efficiency of our method. We will include this result in the Supplementary Information after revision.\n\nFor **parameter sensitivity and selection**, we will include a parameter analysis section in the supplementary information. We analyze the effects of two main parameters: the step size $\\gamma$ and the number of iterations $T$ on the performance of our HeO method. Since revisions are not allowed during the rebuttal phase, we present the results in Figure 2 in the ***PDF file of the global author rebuttal***. This figure shows that HeO performs well across a wide range of step sizes $\\gamma$, provided that the number of iterations $T$ is sufficient. This indicates that HeO is relatively insensitive to the step size $\\gamma$.\n\n\n> 2. The paper demonstrates improved performance over several classic solvers like MDGE, SA, and LQA. However, it does not provide enough evidence or discuss the potential of HeO to achieve state-of-the-art performance on the studied CO problems\n\nOur goal is not to demonstrate that HeO achieves state-of-the-art performance for specific combinatorial optimization tasks, as metaheuristics can be tailored to achieve state-of-the-art results for particular problems [4]. Instead, we aim to highlight the generality of HeO: *it employs a different theoretical approach while performing competitively across a wide range of combinatorial optimization problems without the need for specialized design for each case*. This makes HeO a promising framework for addressing a broad spectrum of combinatorial optimization challenges.\n\nTo achieve this, we have compared our HeO framework against various types of combinatorial optimization problems, including max-cut, 3-SAT, ternary network training, variable selection, and minimum vertex cover. We evaluated HeO alongside advanced methods proposed in recent years, including state-of-the-art higher-order Ising machines [5] and coherent Ising machines (CIMs) and their variants, which are widely used in industry for various combinatorial optimizations [6].\n\n\n\n> 3. The submission seems to be not incomplete as there lack of appendices or supplementary materials that verify the soundness of the theorems presented in Section 3.\n\n\nThank you for pointing this out. We will include supplementary materials with detailed proofs for Theorems 1–3 and additional implementation details. Since revisions to the paper are not permitted during the rebuttal phase, we provide a brief outline of the proof in the *global author rebuttal*.\n\n[1] Liao F Y, Ding L, Zheng Y. Error bounds, PL condition, and quadratic growth for weakly convex functions, and linear convergences of proximal point methods[C]//6th Annual Learning for Dynamics & Control Conference. PMLR, 2024: 993-1005.\n\n[2] Evans L C. Partial differential equations[M]. American Mathematical Society, 2022.\n\n[3] Atenas F, Sagastizábal C, Silva P J S, et al. A unified analysis of descent sequences in weakly convex optimization, including convergence rates for bundle methods[J]. SIAM Journal on Optimization, 2023, 33(1): 89-115.\n\n[4] Bybee C, Kleyko D, Nikonov D E, et al. Efficient optimization with higher-order Ising machines[J]. Nature Communications, 2023, 14(1): 6033\n\n[5] Wang J, Ebler D, Wong K Y M, et al. Bifurcation behaviors shape how continuous physical dynamics solves discrete Ising optimization[J]. Nature Communications, 2023, 14(1): 2510.\n\n[6] Glover, Fred W., and Gary A. Kochenberger, eds. Handbook of metaheuristics. Vol. 57. Springer Science & Business Media, 2003.", "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": "jVNVoRgtB4", "reviewer": "Reviewer_VgKr", "summary": "This work solves combinatorial optimization problems using the gradient method by transforming the discrete problem into a continuous problem. Under the invariant of the optimal solution, the authors transformed the hard continuous problem into an easier problem by changing the objective function using a heating equation and improved the calculation process of the gradient which makes the problem more tractable.", "questions": "1. The discussion of equations (8) to (9) is a little confusing.  Could you please explain how to calculate the projection map in equation(9)?", "limitations": "1. The method does not fully discuss the combinatorial problem with constraints which is the main part of the combinatorial problem.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The authors established the basic theory for the method proposed in the paper.\n2. Extensive experients have been conducted and clearly figures have been presented.", "weaknesses": "1. The motivation of this paper seems to improve the scope of the search-based combinatorial optimization solver, but the proposed method is to change the objective function to improve the efficiency of the gradient method. The method seems a little bit irrelevant to the original motivation of the paper.\n2. The paper does not fully discuss the combinatorial optimization problem with constraints ( only do some experiments on the minimum vertex cover problem ), and the description of the violation of constraints is not sufficient.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 1. The motivation of this paper seems to improve the scope of the search-based combinatorial optimization solver, but the proposed method is to change the objective function to improve the efficiency of the gradient method. The method seems a little bit irrelevant to the original motivation of the paper.\n\nThank you for pointing this out. We will clarify this in the revised version. The primary reason we focus on gradient methods is that, to address the limitations of search-based combinatorial optimization solvers, we reformulated the combinatorial problems into continuous optimization problems. This reformulation allows us to apply heat diffusion to propagate information across the configuration space. Consequently, gradient methods are a natural choice for solving these continuous optimization problems.\n\nAdditionally, gradient methods can be viewed as a specific type of search-based combinatorial optimization solver, where the gradient provides local information that guides the solver such as Gibbs sampling methods in the direction of improving the solution [1].\n\n> 2. The paper does not fully discuss the combinatorial optimization problem with constraints (only do some experiments on the minimum vertex cover problem), and the description of the violation of constraints is not sufficient.\n\nThank you for pointing this out. We will clarify this in the revised version. In Section 4 of the main paper, we present experiments on large-scale instances of the minimum vertex cover problem, where our HeO consistently finds good solutions without violating any constraints. These instances involve millions of nonlinear constraints (see Eq. 18), with edge and vertex counts reaching up to 50 million (see Table 1). Our experiments demonstrate HeO's ability to handle complex, large-scale problems with a massive number of nonlinear constraints. Additionally, HeO can be integrated with other existing techniques for combinatorial optimization problems with constraints, such as Augmented Lagrangian methods, to achieve better performance [2].\n\n\n> 3. The discussion of equations (8) to (9) is a little confusing. Could you please explain how to calculate the projection map in equation(9)?\n\nThank you for pointing this out. We will clarify this in the revised version. The projection of a point $\\mathbf{x} \\in \\mathbb{R}^n$ onto the region $\\mathcal{I}=[0,1]^n$ is calculated for each coordinate $i$ of $x$ as: \n\\begin{align}\n\\mathrm{Proj}_{\\mathcal{I}}[x]_i =\\min(1, \\max(0, x_i)).\n\\end{align}\n\n[1] Grathwohl W, Swersky K, Hashemi M, et al. Oops i took a gradient: Scalable sampling for discrete distributions[C]//International Conference on Machine Learning. PMLR, 2021: 3831-3841.\n\n[2] Birgin E G, Martínez J M. Practical augmented Lagrangian methods for constrained optimization[M]. Society for Industrial and Applied Mathematics, 2014.", "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:
We sincerely thank all the reviewers for their valuable feedback and suggestions. In addition to responding to each reviewer individually, we have included two figures in the ***PDF file*** of this global author rebuttal. Furthermore, to address concerns about theoretical soundness raised by Reviewer uVQ4, we provide a sketch of the proofs of the main theorems presented in Section 3 below.

Since revisions to the paper are not permitted during the rebuttal phase, we will include both the figures and the detailed proofs of the main theorems in the Supplementary Information after the revision.


###  Sketch of the Proof for Theorem 1

We need to show that for any $\tau > 0$, $u(0, \theta)=h(\theta)$ and $u(\tau, \theta)$ have the same global minimas.
It is straightforward to show that the global minima of $u(0, \theta)$ is also a global minima of $u(\tau, \theta)$ for any $\tau > 0$.
The proof for the converse direction relies on the backward uniqueness of the heat equation [1], which asserts that the initial state of a heat equation can be uniquely determined by its state at a time point $\tau$, provided some mild growth conditions, which are satisfied in our paper. 
Denote $u(\tau,\mathbf{x};{\theta})=\mathbb{E}_{p(\mathbf{z})}[f(\mathrm{sgn}({\theta}-(\mathbf{x}+\sqrt{2\tau}\mathbf{z})))]$, where $p(\mathbf{z})$ obeys standard Gaussian distribution. If $\theta^\ast$ is one of the minimas of $h(\theta)$, and $\hat{\theta}$ is one of the minimas of $u(\tau,\theta)$, it can be proved that
\begin{align*}
u(\tau,\mathbf{x};\hat{{\theta}})=u(\tau,\mathbf{x};{{\theta}}^{\ast}),\quad \mathbf{x}\in \mathbb{R}^n
\end{align*}
Using the backward uniqueness of the heat equation, we have
\begin{align*}
    u(0,\mathbf{x};\hat{{\theta}})= u(0,\mathbf{x};{{\theta}}^{\ast}),\quad \mathbf{x}\in\mathbb{R}^n,
\end{align*}
that is
\begin{align*}
    h(\hat{{\theta}}) = h({{\theta}}^{\ast}).
\end{align*}
As a result, $\hat{{\theta}}$ is the one of minimas of $h({\theta})$.


###  Sketch of the Proof for Theorem 2

We first write the $u$ as 

\begin{align*}
    u(\tau,\theta) 
= \mathbb{E}_{p(\mathbf{x},\mathbf{z})} [f(\mathrm{sgn}(\theta+\sqrt{2\tau}\mathbf{z}-\mathbf{x}))],
\end{align*}

where $\mathbf{z} \sim N(\mathbf{0},I)$ is independent of $\mathbf{x}$. Using the property of multi-dimensional Gaussian integral [2], we have

\begin{align*}
    \mathbb{E}_{p(\mathbf{z})}[f(\mathrm{sgn}({\theta}+\sqrt{2\tau}\mathbf{z}-\mathbf{x}))]=f(\tilde{\mathbf{s}}),
\end{align*}

where $\tilde{\mathbf{s}}$ is a random vector determined by $\mathbf{x}$

\begin{align*}
    \tilde{s}_i = \mathrm{erf}(\frac{\theta_i - x_i}{\sqrt{2\tau}}),
\end{align*}

and $\mathrm{erf}(\cdot)$ is the error function. Therefore, we have

\begin{align*}
    u(\tau,{\theta}) = \mathbb{E}_{p(\mathbf{x})}[f(\mathrm{erf}(\frac{{\theta} - \mathbf{x}}{\sqrt{2\tau}}))]
\end{align*}

where $\mathrm{erf}(\cdot)$  is the element-wise error function.

###  Sketch of the Proof for Theorem 3

Define the square loss of ${\theta}$ as $e({\theta}) = (h({\theta})-h({\theta}^{\ast}))^2$
and the error function 
\begin{align*}
    r(\tau,\mathbf{x};{\theta})  =  u(\tau,\mathbf{x};{\theta}) - u(\tau,\mathbf{x};{\theta}^{\ast})
\end{align*}

Define the energy function of the error function $r(\tau,\mathbf{x};{\theta})$ as
\begin{align*}
    E(\tau;{\theta}) = \int_{\mathbb{R}^n}r^2(\tau,\mathbf{x};{\theta})p(\mathbf{x})d\mathbf{x}.
\end{align*}

Use the Harnack's inequality[3] and integration by parts, we have for $0<\tau_1<\tau_2$
\begin{align*}
  E(\tau_1;{\theta}) \leq E(\tau_2;{\theta}) +  \frac{n}{2}\int_{\tau_1}^{\tau_2}\frac{ E(\tau;{\theta})}{\tau} d\tau.
\end{align*} 
Using the Minkowski inequality on the measure $p(\mathbf{x})$, we have
\begin{align*}
      h({\theta})-h({\theta}^{\ast}) \leq 
        \big(\int_{\mathbb{R}^n}  (f(\mathrm{sgn}({\theta}-\mathbf{x}))-u(\tau_1;\mathbf{x};{\theta}))^2 p(\mathbf{x})d\mathbf{x}\big)^{1/2}
     +
     \big(\int_{\mathbb{R}^n}  (f(\mathrm{sgn}({\theta}^{\ast}-\mathbf{x}))-u(\tau_1;\mathbf{x};{\theta}^{\ast}))^2 d\mathbf{x}\big)^{1/2}+E^{1/2}(\tau_1;{\theta}).
\end{align*}
Using the continuity of the heat operator, given $\epsilon>0$, there exists a $\tau_1>0$, such that
\begin{align*}
\begin{aligned}
    \big(\int_{\mathbb{R}^n}  (f(\mathrm{sgn}({\theta}-\mathbf{x}))-u(\tau_1;\mathbf{x};{\theta}))^2 p(\mathbf{x})d\mathbf{x}\big)^{1/2}
     +\big(\int_{\mathbb{R}^n}  (f(\mathrm{sgn}({\theta}^{\ast}-\mathbf{x}))-u(\tau_1;\mathbf{x};{\theta}^{\ast}))^2p(\mathbf{x}) d\mathbf{x}\big)^{1/2}<\epsilon.
\end{aligned}
\end{align*}    
 We then have the error control for $e({\theta})$:
\begin{align*}
    e^{1/2}({\theta}) \leq E^{1/2}(\tau_1;{\theta}) + \epsilon \leq  \big(E(\tau_2;{\theta}) +\frac{n}{2}\int_{\tau_1}^{\tau_2}\frac{ E(\tau;{\theta})}{\tau} d\tau \big)^{1/2}+ \epsilon.
\end{align*}
Noticed that
\begin{align*}
    E(\tau;{\theta})\leq (\breve{f}-f^{\ast})(u(\tau,{\theta}^{\ast})-u(\tau,{\theta})),
\end{align*}
where $ \breve{f}=\max_{\mathbf{s}} f(\mathbf{s}),f^\ast=\min_{\mathbf{s}} f(\mathbf{s})$, and we prove the theorem.

[1] Jie Wu and Liqun Zhang. Backward uniqueness for general parabolic operators in the whole space. Calculus of Variations and Partial Differential Equations, 58:1–19, 2019.

[2] Mobahi H, Fisher J W. On the link between gaussian homotopy continuation and convex envelopes[C]//Energy Minimization Methods in Computer Vision and Pattern Recognition: 10th International Conference, EMMCVPR 2015, Hong Kong, China, January 13-16, 2015. Proceedings 10.

[3] Evans L C. Partial differential equations[M]. American Mathematical Society, 2022.

Meta-review: This paper introduces a novel approach to enhancing combinatorial optimization through heat diffusion, effectively transforming the target function to facilitate better information flow back to the solver. The reviewers unanimously found the idea to be novel and promising. During the discussion, the authors provided clarification and additional theoretical analysis that addressed most concerns.
Decision: Accept (poster)