Combinatorial optimization problems are prevalent in various applications, encompassing circuit design [1], 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, 10], simulated bifurcation [11, 12, 13], coherent Ising machine [14, 15], high-order Ising machine [16], and deep learning techniques [17, 18]. 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 a 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 a thorough evaluation of all these solutions expensive. As a result, current approaches are limited to a narrow receptive field, rendering them blind to distant regions in the solution space and heightening the risk of getting trapped in local minimas or areas with bumpy landscapes. Although methods like large neighborhood search [19], variable neighborhood search [20] and path auxiliary sampling [21] 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.

In this study, we approach the prevalent limitation stated above from a 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 to 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 a key in a dark room, as depicted in Fig. 1. Without light, the person is compelled to rely solely on touching his surrounding space. The tactile provides only localized information, leading to 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 a directional beacon. This would significantly enhance navigational efficiency for finding the key.

Motivated by the above metaphor, we propose a 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] guarantees that the original problem’s optima are unchanged under these transformations. Therefore, information of target functions under different heat diffusion transformations can be cooperatively employed for optimize the original problem (Fig. 1). 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 [24], our research further reveals the potential of heat diffusion, a related thermodynamic phenomenon, in enhancing combinatorial optimization.

We will first formulate the general combinatorial optimization problem and then reframe it in a gradient descent-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\in\mathbb{R}^{n}\mapsto\mathbb{R}$ subjecting to a binary constraints

where $\mathbf{s}$ is binary configuration (bits), and $f(\cdot)$ is the target function. Through the transformation $(\mathbf{s}+1)/2$, our definition aligns with that in [26], where elements of $\mathbf{s}$ take $0$ or $1$. Given the advanced development of gradient-based algorithms, we are interested in converting the discrete optimization problem into a differentiable one, thereby enabling gradient descent. To achieve this purpose, we encode the bits $s_{i}$, $i=1,\ldots,n$ as independent Bernoulli variables $p(s_{i}=\pm 1|\bm{\theta})=0.5\pm(\theta_{i}-0.5)$ with $\theta_{i}\in[0,1]$. In this way, we convert the original combinatorial optimization problem into

where $\mathcal{I}:=[0,1]^{n}$, and $h(\bm{\theta})=\mathbb{E}_{p(\mathbf{s}|\bm{\theta})}[f(\mathbf{s})]$. The minima $\bm{\theta}^{\ast}$ of Eq. (2) is $\bm{\theta}^{\ast}=0.5(\mathrm{sgn}(\mathbf{s}^{\ast})+1)$, given that $\mathbf{s}^{\ast}$ is a minima of the original problem Eq. (1). Here, $\mathrm{sgn}(\cdot)$ is the element-wise sign function. Now Eq. (2) can be solved through gradient descent starting from a given initial $\bm{\theta}_{0}$

where $\gamma$ is the learning rate and $T$ is the iteration number. Unfortunately, this yields a probability distribution $p(\mathbf{s}|\bm{\theta})$ over the configuration space $\{-1,1\}^{n}$, instead of a deterministic binary configuration $\mathbf{s}$ as desired. Although we can manually binarize $\bm{\theta}$ through $\mathbf{B}(\bm{\theta}):=\mathrm{sgn}(\bm{\theta}-0.5)$ to get the binary configuration which maximizes probability $p(\mathbf{s}|\bm{\theta})$, the outcome $f(\mathbf{B}(\bm{\theta}))$ may be much higher than $h(\bm{\theta})$, resulting in significant performance degradation [27]. This suggests that a good gradient-based optimizer should efficiently diminish the uncertainty in the output distribution $p(\mathbf{s}|\bm{\theta})$, which can be measured by its total variance

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 a thermodynamic system, where each parameter $\bm{\theta}$ is referred to as a location and is associated with an initial temperature value $-h(\bm{\theta})$, as shown in Fig. 1. Then the optimization procedure can be described as the process that the solver is walking around the parameter space to find the location $\bm{\theta}^{\ast}$ with the highest temperature (or equivalently, the global minima of $h(\bm{\theta})$). 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 $\bm{\theta}^{\ast}$ 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 $\bm{\theta}^{\ast}$.

We apply our HeO to a variety of NP-hard combinatorial optimization problems to demonstrate its broad applicability. Unless explicitly stated otherwise, we employ the $\tau_{t}$ schedule as $\sqrt{\tau_{t}}=:\sigma_{t}=2(1-t/T)$ for HeO throughout this work. The sensitivity of other parameter settings including the step size $\gamma$ and iterations $T$ 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 is not theoretically necessary, as elaborated in Discussion.

Toy example. We consider the following target function

where $\mathrm{sigmoid}(x)=\tfrac{1}{1+e^{-x}}$, and the elements of the network parameters $\mathbf{a}_{1}\in\mathbb{R}^{n}$, $\mathbf{a}_{2}\in\mathbb{R}^{m}$, $W\in\mathbb{R}^{m,n}$ are uniformly sampled from $[-1,1]$ and fixed during optimizing. According to the universal approximation theory [32], $f(\mathbf{s})$ can approximate any continuous function with sufficiently large $m$, thereby representing a general target function. We compare the performance of our HeO with momentum (Alg. 3) against several representative methods: the conventional gradient-based solver MCGE [28] (with or 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 its uncertainty compared to MCGE.

Quadratic unconstrained binary optimization (QUBO). QUBO is the combinatorial optimization problem with quadratic form target function ($J\in\mathbb{R}^{n\times n}$ is a symmetric matrix with zero diagonals)

This corresponds to the case where $K=2$ in Eq. (10). A well-known class of QUBO is 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 $J$ is 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 [12], bSB [13], dSB [13], CIM [35], and SIM-CIM [15] on max-cut problems in the Biq Mac Library [36]¹¹1https://biqmac.aau.at/biqmaclib.html. 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 HeO is superior to other methods in terms of two metrics.

Polynomial unconstrained binary optimization (PUBO). PUBO is a class of combinatorial optimization problems, in which higher-order interactions between bits $s_{i}$ 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 it as a QUBO problem [38], and the one directly solves PUBO [16]. Quadratization may dramatically increases the dimension of the problem, hence brings heavier computational overhead, while our HeO can be directly used for solving PUBO. A well-known class of PUBO is the Boolean $3$-satisfiability ($3$-SAT) problem [27], which involves determining the satisfiability of a Boolean formula over $n$ Boolean variables $b_{1},\ldots,b_{n}$ where $b_{i}\in\{0,1\}$. The Boolean formula is structured in Conjunctive Normal Form (CNF) consisting of $H$ conjunction ($\wedge$) of clauses, and each clause $h$ is a disjunction ($\vee$) of exactly three literals (either a Boolean variable or its negation). An algorithm of $3$-SAT aims to find the Boolean variables that makes as many as clauses satisfied.

To apply our HeO to the 3-SAT, we encode each Boolean variable $b_{i}$ as $s_{i}$, which is assigned with value $1$ if $b_{i}=1$, otherwise $s_{i}=-1$. For a literal, we define a value $c_{h_{i}}$, which is $-1$ if the literal is the negation of the corresponding Boolean variable, otherwise it is $1$. Then finding the Boolean variables that satisfies as many as clauses is equivalent to minimize the target function

This corresponds to the case where $K=3$ in Eq. (10). 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 ²²2https://www.cs.ubc.ca/~hoos/SATLIB/benchm.html. As shown in Fig. 4, our HeO is superior to other methods in attaining higher quality of solutions and finding more the complete satisfiable solution (solutions that satisfying all clauses). Notably, for the cases of 175-250 variables, our HeO is able to find more complete satisfiable solutions, compared to the 3-order OIM, while the 2-order OIM fails to find any complete solutions [16].

Ternary optimization. Neural networks excel in learning and modeling complex functions, but they also bring about considerable computational demands due to their vast number of parameters. A promising strategy to mitigate this issue is quantization, which converts network parameters into discrete values [39]. However, directly training networks with discrete parameters introduces a 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,0,1}$). Supposed that we have an input-output dataset $\mathcal{D}$ generated by a ground-truth ternary single-layer perceptron $\mathbf{y}=\bm{\Gamma}(\mathbf{v};W_{\mathrm{GT}})=\mathrm{Relu}(W_{\mathrm{GT}}\mathbf{v})$, where $\mathrm{Relu}(x)=\max\{0,x\}$, $W_{\mathrm{GT}}\in\{-1,0,1\}^{m\times n}$ is the ground-truth ternary weight parameter, $\mathbf{v}\in\{-1,0,1\}^{n}$ is the input, and $\mathbf{y}$ is the model output. We aim to find the ternary configuration $W\in\{-1,0,1\}^{m\times n}$ minimizing the loss $\mathrm{MSE}(W,\mathcal{D})=\tfrac{1}{\left|\mathcal{D}\right|}\sum_{(\mathbf{v},\mathbf{y})\in\mathcal{D}}\left\|\bm{\Gamma}(\mathbf{v};W)-\mathbf{y}\right\|^{2}$. 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 $\mathbf{s}$ (see Alg. 5, Appendix for details), and the target function is defined as

As shown in Fig. 5, HeO robustly exceeds MCGE under different dataset size $|\mathcal{D}|$ and output size $m$.

Mixed combinatorial optimization. In high-dimensional linear regression, usually only a small fraction of the variables significantly contribute to prediction. Identifying and selecting a subset of variables with strong predictive power—a process known as variable selection—is crucial, as it 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 a practical alternative, regularization methods like Lasso algorithm are commonly employed [40].

Supposed a dataset is generated from a linear model, in which the relation between input $\mathbf{v}$ and output $y$ is $y=\bm{\beta}^{\ast}\cdot\mathbf{v}+\epsilon$, where $\bm{\beta}^{\ast}$ is the ground-truth linear coefficient and $\epsilon$ is independent Gaussian noise with standard deviation $\sigma_{\epsilon}$. Suppose that only a small proportion (denoted as $q\in(0,1)$) of coordinates in $\beta_{i}^{\ast}$ are non-zero. Our goal is to identify these coordinates through an indicator vector $\mathbf{s}\in\{-1,1\}^{n}$ ($1$ for selection and $-1$ for non-selection) and estimate these non-zero coefficients. The target function of the problem is ($\mathbf{1}\in\mathbb{R}^{n}$ is the all-one vector)

We solve the problem through HeO (Alg. 6, Appendix), where we minimize the loss relative to $\bm{\theta}$ while slowing varying $\bm{\beta}$ via its error gradient. After obtaining the indicator $\mathbf{s}$, we conduct an ordinary least squares regression on the variables selected by the $\mathbf{s}$ to estimate their coefficients, and treat other variables’ coefficients as zero. As shown in Fig.6, our HeO outperforms both Lasso regression and the more advanced $L_{0.5}$ regression [41] in terms of producing more accurate indicators $\mathbf{s}$ and achieving lower test prediction errors across various $q$ and $\sigma_{\epsilon}$ settings. Importantly, to give the variable selection prediction, our HeO does not need to know the level of $q$ and $\sigma_{\epsilon}$ in advance.

Constrained binary optimization.

Minimum vertex cover (MVC) is a class of the constrained combinatorial optimization which has wide applications [44], as illustrated in Fig. 7. Given an undirected graph $\mathcal{G}$ with vertex set $\mathcal{V}$ and edge set $\mathcal{E}$, the MVC is to find the minimum subset $\mathcal{V}_{c}\subset\mathcal{V}$, so that for each edge $e\in\mathcal{E}$, at least one of its endpoints belongs to $\mathcal{V}_{c}$. The target function and the constrains are expressed as

We combine the HeO with penalty function (Alg. 7, Appendix) for solving MVC. We compare our HeO with FastVC [45], a powerful MVC heuristic algorithm, on massive real world graph datasets ³³3http://networkrepository.com/. For a fair comparison, we keep the run time of two algorithms as the same for each dataset. As shown in Tab. 1, our HeO can find smaller cover sets than that of FastVC.

Existing model-based combinatorial optimization approaches encode the solution space via a parameterized 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-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 distribution. In contrast, HeO needs no training, and it is unnecessary to strictly adhering to the reverse time sequence $\tau_{t}$ in the optimization process, as under different $\tau$, the function $u(\tau,\bm{\theta})$ shares the same optima with that of the original problem $h(\bm{\theta})$. This claim is empirically corroborated in Fig. S3, Appendix, where HeO applying non-monotonic schedules of $\tau_{t}$ still demonstrates superior performance.

Our HeO can be viewed as a stochastic Gaussian continuation (GC) method [49] with projection, which 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(\bm{\theta},\tau_{t})$) at each $t$ up to some criteria, whereas HeO merely performs a single-step gradient descent. Also, Eq. (8) corresponds to a variation of the evolution strategy (ES), a robust optimizer for non-differentiable function [51], while HeO use a different gradient estimation (Eq. (11)), 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 $\tau$, the smoothed function $u(\tau,\bm{\theta})$ retains the optima of the original function $h(\bm{\theta})$ (Thm. 1). This distinguishes our HeO from methods based on quantum adiabatic theory [9], bifurcation theory [12] and other relaxation strategies [54], in which the optima of the smoothed function can be different from the original one [27, 55]. This is empirically verified in Fig. S3, Appendix.

The heat equation in our HeO can be naturally extended to general parabolic differential equations, given that a broad spectrum of them obey the backward uniqueness [56]. For example, we can use $\partial_{\tau}u(\tau,\bm{\theta})=\bigtriangledown_{\bm{\theta}}[A\bigtriangledown_{\bm{\theta}}u(\tau,\bm{\theta})]$ to replace the Eq. (6), where $A$ is a 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 $\tau_{t}$ does not necessarily have to be monotonic due to the cooperative optimization property of HeO, it is feasible to explore various $\tau_{t}$ schedules (possibly non-monotonic) to further enhance performance. Moreover, our HeO can be integrated with other existing techniques for combinatorial optimization problems with constraints, such as 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 it is cumbersome to encode integer variables through the Bernoulli distribution in our framework. Nevertheless, integrating HeO with other techniques such as advanced Metropolis-Hastings algorithm [60] may path to broaden the applicability of our methodology to a wider range of combinatorial optimization problems. Besides, our HeO allows for further customization by incorporating additional terms that integrate problem-specific prior knowledge or 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]:

where 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 [62].

In conclusion, grounded in the heat diffusion, we present a framework called Heat diffusion optimization (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.