Efficient Combinatorial Optimization via Heat Diffusion

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.

Thank 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.

Additionally, 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].

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.

Thank 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].

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)?

Thank 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:

$

\mathrm{Proj}_{\mathcal{I}}[x]_i =\min(1, \max(0, x_i)).

$

[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.

[2] Birgin E G, Martínez J M. Practical augmented Lagrangian methods for constrained optimization[M]. Society for Industrial and Applied Mathematics, 2014.

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.

Thank you for the valuable suggestions.

For 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]:

$

\left| \bigtriangledown_{{\theta}}u(\tau,{\theta} )\right| \leq \frac{C}{\sqrt{\tau}},

$

where $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].

For 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.

For 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$.

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

Our 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.

To 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].

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.

Thank 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.

[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.

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

[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.

[4] Bybee C, Kleyko D, Nikonov D E, et al. Efficient optimization with higher-order Ising machines[J]. Nature Communications, 2023, 14(1): 6033

[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.

[6] Glover, Fred W., and Gary A. Kochenberger, eds. Handbook of metaheuristics. Vol. 57. Springer Science & Business Media, 2003.

Thank you for your constructive comments. There is considerable flexibility and feasibility in customizing HeO for specific problems.

First, as HeO is a gradient-based optimizer, it can be tailored to specific problems by designing more refined step schedules, such as adaptive step rules. Additionally, leveraging second-order optimization techniques like momentum (which we utilized in our paper) or Adam, as well as ensemble methods, can improve the overall performance of HeO on particular problem instances.

Second, as discussed in Line 265, Section 5, we can customize HeO by designing a preconditioned matrix $A$ to reshape the heat equation. Prior knowledge about the problem can be embedded within the structure of $A$, such as by accounting for the relative importance between different dimensions of the discrete configuration $\mathbf{s}$ or by setting $A$ based on the Fisher information matrix of the parameter $\theta$. This approach can lead to a natural gradient descent method, enhancing the efficiency of the optimization process.

Third, HeO allows for further customization by integrating problem-specific prior knowledge directly into the target function. By adding extra terms, we can improve the loss landscape or guide the search direction to meet particular purposes, thereby improving the quality of the solutions found.

Fourth, HeO can be hybridized with other metaheuristic algorithms to explore the configuration space more effectively. Specifically, we can iteratively refine the solution by alternating between HeO and other metaheuristic update rules.

We will incorporate these discussions into Section 5 in the revised version of the paper.

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.

Good 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.

[1] Sun H, Dai H, Xia W, et al. Path auxiliary proposal for MCMC in discrete space[C]//International Conference on Learning Representations. 2021.

1. Could you elaborate on how K (Page 4, Line 129) is selected in different scenarios?

Thank 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:

For minimum vertex cover problem, where $f(\mathbf{s})$ is a linear function (Eq. 18), we have K=1.

For Quadratic unconstrained binary optimization, where $f(\mathbf{s})$ is a quadratic function (Eq. 14), we have K=2.

For 3-satisfiability problem, where $f(\mathbf{s})$ is a cubic function (Eq. 15), we have K=3.

2. Why does the energy in Fig. 2 not monotonically decrease?

In simulated annealing, non-monotonicity arises due to the stochastic nature of the annealing process, which occasionally permits transitions from lower to higher energy states.

In 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).

3. The algorithmic complexity, computational time, and the scalability of the HeO, especially in very large-dimensional spaces.

Thank 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.

4. Is it possible to apply your heat diffusion framework to other types of combinatorial optimization problems?

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

5. More baselines are needed.

We 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.

6. On Page 3, line 105, the definitions of $\Delta_\theta$ and the function $u$ are not given.

Thank 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}$.

7. The ability of global search is not theoretically analyzed.

Thank 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]:

$

\left| \bigtriangledown_{{\theta}}u(\tau,{\theta} )\right| \leq \frac{C}{\sqrt{\tau}},

$

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

[1] Korte B H, Vygen J, Korte B, et al. Combinatorial optimization[M]. Berlin: Springer, 2011.

[2] Bybee C, Kleyko D, Nikonov D E, et al. Efficient optimization with higher-order Ising machines[J]. Nature Communications, 2023, 14(1): 6033.

[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.

[4] Glover, Fred W., and Gary A. Kochenberger, eds. Handbook of metaheuristics. Vol. 57. Springer Science & Business Media, 2003.

[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.

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

[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.