Controlling Continuous Relaxation for Combinatorial Optimization

Authors: yuma  ichikawa
fujitsu  limited
japan  kanagawa
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 12:15:47

Abstract: Unsupervised learning (UL)-based solvers for combinatorial optimization (CO) train a neural network that generates a soft solution by directly optimizing the CO objective using a continuous relaxation strategy. These solvers offer several advantages over traditional methods and other learning-based methods, particularly for large-scale CO problems. However, UL-based solvers face two practical issues: (I) an optimization issue, where UL-based solvers are easily trapped at local optima, and (II) a rounding issue, where UL-based solvers require artificial post-learning rounding from the continuous space back to the original discrete space, undermining the robustness of the results. This study proposes a Continuous Relaxation Annealing (CRA) strategy, an effective rounding-free learning method for UL-based solvers. CRA introduces a penalty term that dynamically shifts from prioritizing continuous solutions, effectively smoothing the non-convexity of the objective function, to enforcing discreteness, eliminating artificial rounding. Experimental results demonstrate that CRA significantly enhances the performance of UL-based solvers, outperforming existing UL-based solvers and greedy algorithms in complex CO problems. Additionally, CRA effectively eliminates artificial rounding and accelerates the learning process.

Main paper: 
 introduction: 
 The objective of combinatorial optimization (CO) problems is to find the optimal solution from a discrete space, and these problems are fundamental in many real-world applications [Papadimitriou and Steiglitz, 1998] . Most CO problems are NP-hard or NP-complete; making it challenging to solve large-scale problems within feasible computational time. Traditional methods frequently depend on heuristics to find approximate solutions; however, considerable insights into the specific problems are required. Alternatively, problems can be formulated as integer linear programming (ILP) and solved using ILP solvers. However, ILP lacks scalability when applied to problems with graph structures.Recently, several studies have used machine learning methods to handle CO problems by learning heuristics. Most of these studies focus on supervised learning (SL)-based solvers [Hudson et al., 2021 , Joshi et al., 2019 , Gasse et al., 2019 , Selsam et al., 2018 , Khalil et al., 2016] , which require optimal solutions to CO problems as supervision during training. However, obtaining optimal solutions is challenging in practice, and SL-based solvers often fail to generalize well [Yehuda et al., 2020] . Reinforcement learning (RL)-based solvers [Yao et al., 2019 , Chen and Tian, 2019 , Yolcu and Póczos, 2019 , Nazari et al., 2018 , Khalil et al., 2017 , Bello et al., 2016] avoid the need for optimal solutions but often suffer from notoriously unstable training due to poor gradient estimation and hard explorations [Mnih et al., 2015 , Tang et al., 2017 , Espeholt et al., 2018] . Unsupervised learning (UL)-based solvers have recently attracted much attention [Schuetz et al., 2022a , Karalias and Loukas, 2020 , Amizadeh et al., 2018] . UL-based solvers follow a continuous relaxation approach, training a UL model to output a soft solution to the relaxed CO problem by directly optimizing a differentiable objective function, offering significantly stable and fast training even for large-scale CO problems. Notably, the physics-inspired GNN (PI-GNN) solver [Schuetz et al., 2022a] employs graph neural networks (GNN) to automatically learn instance-specific heuristics and performs on par with or outperforms existing solvers for CO problems with millions of variables without optimal solutions. While these offer some advantages over traditional and other machine learning-based solvers, they face two practical issues. The first issue is an optimization issues where UL-based solvers are easily trapped at local optima. Due to this issue, Angelini and Ricci-Tersenghi [2023] demonstrated that the PI-GNN solver [Schuetz et al., 2022a] could not achieve results comparable to those of the degree-based greedy algorithm (DGA) [Angelini and Ricci-Tersenghi, 2019] on maximum independent set (MIS) problems in random regular graphs (RRG) . Wang and Li [2023] also pointed out the importance of using dataset or history, and initializing the GNN with outputs from greedy solvers to help the PI-GNN solver overcome optimization challenges. This issue is a crucial bottleneck to the applicability of this method across various real-world applications. The second issue relates to the inherent ambiguity of the continuous relaxation approach. This approach necessitates artificial rounding from the soft solution, which may include continuous values, back to the original discrete solution, potentially undermining the robustness of the results. While linear relaxation can provide an optimal solution for original discrete problems on bipartite graphs [Hoffman and Kruskal, 2010] , it typically leads to solutions with 1/2 values, which is known to half-integrality [Nemhauser and Trotter Jr, 1974] , in which existing rounding methods [Schuetz et al., 2022b , Wang et al., 2022] completely lose their robustness. For NP-hard problems with graph structures, such as the MIS and MaxCut, semidefinite programming (SDP) relaxations have been proposed as effective approximation methods [Lovász, 1979, Goemans and Williamson, 1995] . However, these approaches rely on rounding techniques, such as spectral clustering [Von Luxburg, 2007] , to transform relaxed solutions into feasible ones, which often fails to obtain optimal solutions.To address these issues, we propose the Continuous Relaxation Annealing (CRA). CRA introduces a penalty term to control the continuity and discreteness of the relaxed variables, with a parameter γ to regulate the intensity of this penalty term. When the parameter γ is small, the relaxed variable tends to favor continuous solutions, whereas a large γ biases them toward discrete values. This penalty term also effectively eliminates local optimum. Moreover, a small γ forces the loss function to approach a simple convex function, encouraging active exploration within the continuous space. CRA also includes an annealing process, where γ is gradually increased until the relaxed variables approach discrete values, eliminating the artificial rounding from the continuous to the original discrete space after learning. In this study, the solver that applies the CRA to the PI-GNN solver is referred to as the CRA-PI-GNN solver. We also demonstrate the benefits of the CRA through experiments on benchmark CO problems, including MIS, maximum cut (MaxCut), and diverse bipartite matching (DBM) problems across graphs of varying sizes and degrees. The experimental results show that the CRA significantly enhances the performance of the PI-GNN solver, outperforming the original PI-GNN solver, other state-of-the-art learning-based baselines, and greedy algorithms. This improvement is achieved by directly optimizing each instance without any history, e.g., previous optimal solutions and the information of other instances. Additionally, these experiments indicate that the CRA accelerates the learning process of the PI-GNN solver. Notably, these results overcome the limitations pointed out by Angelini and Ricci-Tersenghi [2023] , Wang and Li [2023] , highlighting the further potential of UL-based solvers.Notation We use the shorthand expression [N ] = {1, 2, . . . , N }, where N ∈ N. I N ∈ R N ×N denotes an N × N identity matrix, 1 N denotes the vector (1, . . . , 1) ⊤ ∈ R N , and 0 N denotes the vector (0, . . . , 0) ⊤ ∈ R N . G(V, E) represents an undirected graph, where V is the set of nodes with cardinality |V | = N , and E ⊆ V × V denotes the set of edges. For a graph G(V, E), A ij denotes the adjacency matrix, where A ij = 0 if an edge (i, j) does not exist and A ij > 0 if the edge is present. 
background: 
 Combinatorial optimization The goal of this study is to solve the following CO problem.EQUATIONwhere C ∈ C denotes instance-specific parameters, such as a graph G = (V, E), and C represents the set of all possible instances. f : X × C → R denotes the cost function. Additionally, x = (x i ) 1≤i≤N ∈ {0, 1} N is a binary vector to be optimized, and X (C) ⊆ {0, 1} N denotes the feasible solution space, typically defined by the following equality and inequality constraints.X (C) = {x ∈ {0, 1} N | ∀i ∈ [I], g i (x; C) ≤ 0, ∀j ∈ [J], h j (x; C) = 0}, I, J ∈ N,Here, for i ∈ [I], g i : {0, 1} N × C → R denotes the inequality constraint, and for j ∈ [J], h j : {0, 1} N × C → R denotes the equality constraint. Following UL-based solvers [Wang et al., 2022 , Schuetz et al., 2022a , Karalias and Loukas, 2020] , we reformulate the constrained problem into an equivalent unconstrained form using the penalty method [Smith et al., 1997] :EQUATIONHere, for all i ∈ [I + J], v : {0, 1} N × C → R is the penalty term, which increases when the constraints are violated. For example, the penalty term is defined as follows:∀i∈ [I], j ∈ [J], v i (x; C) = max(0, g i (x; C)), ∀j ∈ [J], v j (x; C) = (h j (x; C)) 2and λ = (λ i ) 1≤i≤I+J ∈ R I+J denotes the penalty parameters that control the trade-off between constraint satisfaction and cost optimization. Note that, as λ increases, the penalty for constraint violations becomes more significant. In the following, we provide an example of this formulation.The MIS problem is a fundamental NP-hard problem [Karp, 2010] , defined as follows. Given an undirected graph G(V, E), an independent set (IS) is a subset of nodes I ∈ V where any two nodes are not adjacent. The MIS problem aims to find the largest IS, denoted as I * . In this study, ρ denotes the IS density, defined as ρ = |I|/|V |. Following Schuetz et al.[2022a], a binary variable x i is assigned to each node i ∈ V . The MIS problem can be formulated as follows:EQUATIONwhere the first term maximizes the number of nodes assigned a value of 1, and the second term penalizes adjacent nodes assigned 1 according to the penalty parameter λ.Learning for CO problems involves training an algorithm A θ (•) : C → {0, 1} N parameterized by a neural network (NN), where θ denotes the parameters. For a given instance C ∈ C, this algorithm generates a valid solution x = A θ (C) ∈ X (C) and aims to minimize f ( x; C). Several approaches have been proposed to train A θ . This study focuses on UL-based solvers, which do not use a labeled solution x * ∈ argmin x∈X (C) f (x; C) during training [Wang et al., 2022 , Schuetz et al., 2022a , Karalias and Loukas, 2020 , Amizadeh et al., 2018] . In the following, we outline the details of the UL-based solvers.The UL-based solvers employ a continuous relaxation strategy to train NN. This continuous relaxation strategy reformulates a CO problem into a continuous optimization problem by converting discrete variables into continuous ones. A typical example of continuous relaxation is expressed as follows:min p l(p; C, λ), l(p; C, λ) = ∆ f (p; C) + m+p i=1 λ i vi (p; C),where p = (p i ) 1≤i≤N ∈ [0, 1] N represents a set of relaxed continuous variables, where each binary variable Wang et al. [2022] and Schuetz et al. [2022a] formulated A θ (C) as the relaxed continuous variables, defined as A θ (•) : C → [0, 1] n . In the following discussions, we denote A θ as p θ to make the parametrization of the relaxed variables explicit. Then, p θ is optimized by directly minimizing the following label-independent function:x i ∈ {0, 1} is relaxed to a continuous counterpart p i ∈ [0, 1], and f : [0, 1] N × C → R denotes the relaxed form of f such that f (x; C) = f (x; C) for x ∈ {0, 1} N . The relation between each constraint v i and its relaxation vi is similar for i ∈ [I +J], meaning that ∀i ∈ [I +J], vi (x; C) = v i (x; C) for x ∈ {0, 1} N .EQUATIONAfter training, the relaxed solution p θ is converted into discrete variables using artificial rounding p θ , where ∀i ∈ [N ], x i = int(p θ,i (C)) based on a threshold [Schuetz et al., 2022a] , or alternatively, a greedy method [Wang et al., 2022] . Two types of schemes for UL-based solvers have been developed based on this formulation.(Type I) Learning generalized heuristics from history/data One approach, proposed by Karalias and Loukas [2020] , aims to automatically learn effective heuristics from historical dataset instances D = {C µ } P µ=1 and then apply these learned heuristics to a new instance C * , through inference. Note that this method assumes that either the training dataset is easily obtainable or that meaningful data augmentation is feasible. Specifically, given a set of training instances D = (C µ ), sampled independently and identically from a distribution P (C), the goal is to minimize the average loss function min θ P µ=1 l(θ; C µ , λ). However, this method does not guarantee quality for a test instance, C * . Even if the training instances D are extensive and the test instance C follows P (C), low average performance E C∼P (C) [ l(θ; C)] may not guarantee a low l(θ; C) for on a specific C. To address this issue, Wang and Li [2023] introduced a meta-learning approach where NNs aim to provide good initialization for future instances rather than direct solutions.(Type II) Learning effective heuristics on a specific single instance Another approach, known as the PI-GNN solver [Schuetz et al., 2022a,b] , automatically learns instance-specific heuristics for a single instance using the instance parameter C by directly applying Eq. ( 4). This approach addresses CO problems on graphs, where C = G(V, E), and employs GNNs for the relaxed variables p θ (G).Here, an L-layered GNN is trained to directly minimize l(θ; C, λ), taking as input a graph G and the embedding vectors on its nodes, and outputting the relaxed solution p θ (G) ∈ [0, 1] N . A detailed description of GNNs is provided in Appendix E.2. Note that this setting is applicable even when the training dataset D is difficult to obtain. The overparameterization of relaxed variables is expected to smooth the objective function by introducing additional parameters to the optimization problem, similar to the kernel method. However, minimizing Eq. 4 for a single instance can be time-consuming compared to the inference process. Nonetheless, for large-scale CO problems, this approach has been reported to outperform other solvers in terms of both computational time and solution quality [Schuetz et al., 2022a,b] .Note that, while both UL-based solvers for multiple instances (Type I) and individual instances (Type II) are valuable, this study focuses on advancing the latter: a UL-based solver for a single instance. Both types of solvers are applicable to cost functions that meet a particular requirement due to their reliance on a gradient-based algorithm to minimize Eq (4).Assumption 2.1 (Differentiable cost function). The relaxed loss function l(θ; C, λ) and its partial derivative ∂ l(θ;C,λ) /∂θ are accessible during the optimization process.These requirements encompass a nonlinear cost function and interactions involving many-body interactions, extending beyond simple two-body interactions. 
continuous relaxation annealing for ul-based solvers: 
 In this section, we discuss the practical issues associated with UL-based solvers and then introduce continuous relaxation annealing (CRA) as a proposed solution.UL-based solvers (Type II) [Schuetz et al., 2022a,b] are effective in addressing large-scale CO problems. However, these solvers present following two practical issues, highlighted in several recent studies [Wang and Li, 2023, Angelini and Ricci-Tersenghi, 2023] . Additionally, we numerically validate these issues; see Appendix F.1 for detailed results.(I) Ambiguity in rounding method after learning UL-based solvers employ a continuous relaxation strategy to train NNs and then convert the relaxed continuous variables into discrete binary values through artificial rounding as discussed in Section 2.1. This inherent ambiguity in continuous relaxation strategy often results in potential discrepancies between the optimal solutions of the original discrete CO problem and those of the relaxed continuous one. Continuous relaxation expands the solution space, often producing continuous values that lower the cost compared to an optimal binary value. Indeed, while linear relaxation can provide an optimal solution for discrete problems on bipartite graphs [Hoffman and Kruskal, 2010] , it typically results in solutions with 1 /2 values, which is known to half-integrality [Nemhauser and Trotter Jr, 1974] . Existing rounding methods [Schuetz et al., 2022b , Wang et al., 2022] often lose robustness in these scenarios. In practice, PI-GNN solver often outputs values near 1 /2, underscoring the limitations of current rounding techniques for UL-based solvers.(II) Difficulty in optimizing NNs Recently, Angelini and Ricci-Tersenghi [2023] demonstrated that PI-GNN solver falls short of achieving results comparable to those of the degree-based greedy algorithm (DGA) [Angelini and Ricci-Tersenghi, 2019] when solving the MIS problems on RRGs. Angelini and Ricci-Tersenghi [2023] further emphasized the importance of evaluating UL-based solvers on complex CO problems, where greedy algorithms typically perform worse. A representative example is the MIS problems on RRGs with a constant degree d > 16, where a clustering transition in the solution space creates barriers that impede optimization. Moreover, Wang and Li [2023] emphasized the importance of using training/historical datasets, D = {C µ } 1≤µ≤P , which contain various graphs and initialization using outputs from greedy solvers, such as DGA and RGA for MIS problems. Their numerical analysis indicated that PI-GNN solver tends to get trapped in local optima when directly optimized directly for a single instance without leveraging a training dataset D. However, in a practical setting, systematic methods for generating or collecting training datasets D to effectively avoid local optima remains unclear. Additionally, training on instances that do not contribute to escaping local optima is time-consuming. Therefore, it is crucial to develop an effective UL-based solver that can operate on a single instance without relying on training data, D. Our numerical experiments, detailed in Appendix F.1, also confirmed this optimization issue. They demonstrated that as problem complexity increases, the PI-GNN solver is often drawn into trivial local optima, p θ = 0 N , in certain problems. This entrapment results in prolonged plateaus that significantly slow down the learning process and, in especially challenging cases, can render learning entirely infeasible. Our numerical experiments, detailed in Appendix F.1, also validated this optimization issue, demonstrating that as the problem complexity increases, PI-GNN solver tends to be absorbed into the trivial local optima p θ = 0 N in some problems, resulting in prolonged plateaus which significantly decelerates the learning process and, in particularly challenging cases, can render learning entirely infeasible.Penalty term to control discreteness and continuity To address these issues, we propose a penalty term to control the balance between discreteness and continuity in the relaxed variables, formulated as follows:r(p; C, λ, γ) = l(p; C, λ) + γΦ(p), Φ(p) = ∆ N i=1 (1 -(2p i -1) α ), α ∈ {2n | n ∈ N + }, (5)where γ ∈ R is a penalty parameter, and the even number α denote a curve rate. When γ is negative, i.e., γ < 0, the relaxed variables tend to favor the continuous space, smoothing the non-convexity of the objective function l(p; C, λ) due to the convexity of the penalty term Φ(p). In contrast, when γ is positive, i.e., γ > 0, the relaxed variables tend to favor discrete space, smoothing out the continuous solution into discrete solution. Formally, the following theorem holds as λ approaches ±∞.Theorem 3.1. Assuming the objective function l(p; C) is bounded within the domain [0, 1] N , as γ → +∞, the relaxed solutions p * ∈ argmin p r(p; C, λ, γ) converge to the original solutions x * ∈ argmin x l(x; C, λ). Moreover, as γ → -∞, the loss function r(p; C, λ, γ) becomes convex, and the relaxed solution 1 N /2 = argmin p r(p, C, λ, γ) is unique.For the detailed proof, refer to Appendix B.1. Theorem 3.1 can be generalized for any convex function Φ(p; C) that has a unique maximum at 1 N /2 and achieves a global minimum for all p ∈ {0, 1} N ; an example is binary cross entropy Sanokowski et al. [2024] for the UL-based solvers (Type I). Additionally, the penalty term eliminates the stationary point p * = 0 N described in Section 3.1, preventing convergence to a plateau. For UL-based solvers, the penalty term is expressed as follows: As γ increases, it promotes optimal discrete solutions by smoothing away suboptimal continuous ones.Φ CE (p) = N i=1 (p i log p i +(1-p i ) log(1-p i )), introduced by Sun et al. [2022],EQUATIONWe propose an annealing strategy that gradually anneals the penalty parameter γ in Eq. ( 6). Initially, a negative gamma value, i.e., γ < 0, is chosen to leverage the properties, facilitating broad exploration by smoothing the non-convexity of l(θ; C, λ) and eliminating the stationary point p * = 0 N to avoid the plateau, as discussed in Section 3.1. Subsequently, the penalty parameter γ is gradually increased to a positive value, γ > 0, with each update of the trainable parameters (one epoch), until the penalty term approaches zero, i.e., Φ(θ, C) ≈ 0, to automatically round the relaxed variables by smoothing out suboptimal continuous solutions oscillating between 1 or 0. A conceptual diagram of this annealing process is shown in Fig. 1 .Note that employing the binary cross-entropy Φ CE (p) is infeasible for UL-based solvers when γ > 0, as the gradient ∂ΦCE(p) /∂pi diverges to ±∞ at 0 or 1. In deed, when γ = 0, most relaxed variables typically approach binary values, with a relatively small number of variables oscillating between 0 and 1. This gradient divergence issue in Φ CE (p) makes the learning infeasible without additional techniques, such as gradient clipping. In contrast, the gradient of the penalty term in Eq. 5, ∂Φ(p) /∂pi, is bounded within [-2α, 2α] for any γ, preventing the gradient divergence issue seen in Φ CE (p). Additionally, by increasing α, the absolute value of the gradient near 1 /2 becomes smaller, allowing for control over the smoothing strength toward a discrete solution near 1 /2.We also propose an early stopping strategy that monitors both the loss function and the penalty term, halting the annealing and learning processes when the penalty term approaches zero, i.e., Φ(θ; C) ≈ 0. Various annealing schedules can be considered; in this study, we employ the following scheduling: γ(τ + 1) ← γ(τ ) + ε, where the scheduling rate ε ∈ R + is a small constant, and τ denotes the update iterations of the trainable parameters. We refer to the PI-GNN solver with this continuous relaxation annealing as CRA-PI-GNN solver. Here, two additional hyperparameters are introduced: the initial scheduling value γ(0) and the scheduling rate ε. Numerical experiments suggest that better solutions are obtained when γ(0) is set to a small negative value and ε is kept low. The ablation study are presented in Appendix F.5. 
related work: 
 Previous works on UL-based solvers have addressed various problems, such as MaxCut problems [Yao et al., 2019] and traveling salesman problems [Hudson et al., 2021] , using carefully tailored problemspecific objectives. Some studies have also explored constraint satisfaction problems [Amizadeh et al., 2018 , Toenshoff et al., 2019] , but applying these approaches to broader CO problems often requires problem-specific reductions. Karalias and Loukas [2020] proposed Erdős Goes Neural (EGN) solver, an UL-based solver for general CO problems based on Erdős' probabilistic method. This solver generate solutions through an inference process using training instances. Subsequently, Wang et al. [2022] proposed an entry-wise concave continuous relaxation, broadening the EGN solver to a wide range of CO problems. In contrast, Schuetz et al. [2022a,b] proposed PI-GNN solver, an UL-based solver for a single CO problems that automatically learns problem-specific heuristics during the training process. However, Angelini and Ricci-Tersenghi [2023], Boettcher [2023] pointed out the optimization difficulties where PI-GNN solver failed to achieve results comparable to those of greedy algorithms. Wang and Li [2023] also claimed optimization issues with PI-GNN solver, emphasizing the importance of learning from training data and history to overcome local optima. They then proposed Meta-EGN solvers, a meta-learning approach that updates NN network parameters for individual CO problem instances. Furthermore, to address these optimization issue, Lin et al. [2023] , Sun et al. [2022] , Sanokowski et al. [2024] proposed annealing strategy similar to simulated annealing [Kirkpatrick et al., 1983] . 
experiments: 
 We begin by evaluating the performance of CRA-PI-GNN solver on the MIS and the MaxCut benchmark problems across multiple graphs of varying sizes, demonstrating that CRA effectively overcomes optimization challenges without relying on data/history D. We then extend the evaluation to the DBM problems, showing the applicability to more practical CO problems. For the objective functions and the detailed explanations, refer to Appendix E.1.Baseline methods In all experiments, the baseline methods include the PI-GNN solver [Schuetz et al., 2022a] as the direct baseline of a UL-based solver for a single instance. For the MIS problems, we also consider the random greedy algorithm (RGA) and DGA [Angelini and Ricci-Tersenghi, 2019 ] as heuristic baselines. For the MaxCut problems, RUN-CSP solver [Toenshoff et al., 2019] is considered as an additional baseline, and a standard greedy algorithm and SDP based approximation algorithm [Goemans and Williamson, 1995] are considered as an additional classical baseline. The parameters for the Goemans-Williamson (GW) approximation are all set according to the settings in Schuetz et al. [2022b] . The implementation used the open-source CVXOPT solver with CVXPY [Mehta, 2019] as the modeling interface. Note that we do not consider UL-based solvers for learning generalized heuristics [Karalias and Loukas, 2020 , Wang et al., 2022 , Wang and Li, 2023] , which rely on training instances D = {C µ } P µ=1 . The primary objective of this study is to evaluate whether CRA-PI-GNN solver can surpass the performance of both PI-GNN solver and greedy algorithms. However, for the MIS problem, EGN solver [Karalias and Loukas, 2020] and Meta-EGN solver [Wang and Li, 2023] are considered to confirm that CRA can overcome the optimization issues without training instances.Implementation The objective of the numerical experiments is to compare the CRA-PI-GNN solver with the PI-GNN solver. Thus, we follow the same experimental configuration described as the experiments in Schuetz et al. [2022a] , employing a simple two-layer GCV and GraphSAGE [Hamilton et al., 2017] implemented by the Deep Graph Library [Wang et al., 2019] ; Refer to Appendix D.1 for the detailed architectures. We use the AdamW [Kingma and Ba, 2014] optimizer with a learning rate of η = 10 -4 and weight decay of 10 -2 . The GNNs are trained for up to 5 × 10 4 epochs with early stopping, which monitors the summarized loss function S s=1 l(P :,s ) and the entropy term Φ(P ; γ, α) with tolerance 10 -5 and patience 10 3 epochs; Further details are provided in Appendix D.2. We set the initial scheduling value to γ(0) = -20 for the MIS and matching problems, and we set γ(0) = -6 for the MaxCut problems with the scheduling rate ε = 10 -3 and curve rate α = 2 in Eq. ( 6). These values are not necessarily optimal, and refining these parameters can lead to better solutions; Refer to Appendix F.5 and Appendix F.6 for an ablation study of these parameters. Once the training process is complete, we apply projection heuristics to map the obtained soft solutions back to discrete solutions using simple projection, where for all i ∈ [N ], we map p θ,i into 0 if p θ,i ≤ 0.5 and p θ,i into 1 if p θ,i > 0.5. However, due to the early stopping in Section 3.2, the CRA-PI-GNN solver ensures that for all benchmark CO problems, the soft solution at the end of the training process became 0 or 1 within the 32-bit Floating Point range in Pytorch GPU; thus, it is robust against a given threshold, which we set to 0.5 in our experiments. Additionally, no violations of the constraints were observed in our numerical experiments. Thus, following results presented in are feasible solutions.Evaluation metrics Following Karalias and Loukas [2020] , We use the approximation ratio (ApR), formulated as ApR = f (x; C)/f (x * ; C), where x * is optimal solution. For the MIS problems, we The dashed lines represents the theoretical upper bounds [Parisi, 1980] . evaluate the ApRs using the theoretical optimal cost [Barbier et al., 2013] and the independent set density ρ relative to the theoretical results. For the MaxCut problems on RRGs, we adopt the cut ratio ν against the theoretical upper bound [Parisi, 1980 , Dembo et al., 2017] ; see Appendix E.1 for the details. All the results for the MIS and MaxCut problems are summarized based on 5 RRGs with different random seeds. In the case of the MaxCut Gset problem, the ApR is calculated compared to the known best cost functions. Regarding the DBM problems, we calculate the ApR against the global optimal, identified using Gurobi 10.0.1 solver with default settings.Degree dependency of solutions using CRA First, we compare the performance of the PI-GNN and CRA-PI-GNN solvers using GCV, as in Schuetz et al. [2022a] . Fig. 2 shows the independent set density ρ d as a function of degree d obtained by the PI-GNN and CRA-PI-GNN solvers compared with the theoretical results [Barbier et al., 2013] . Across all degrees d, the CRA-PI-GNN solver outperformed the PI-GNN solver and approached the theoretical results, whereas the PI-GNN solver fail to find valid solutions, especially for d ≥ 15, as pointed by the previous studies [Angelini and Ricci-Tersenghi, 2023, Wang and Li, 2023] .Response to Angelini and Ricci-Tersenghi [2023] and Wang and Li [2023] MIS problems on RRGs with a degree d larger than 16 is known to be hard problems [Barbier et al., 2013] . As discussed in Section 3.1, Angelini and Ricci-Tersenghi [2023] , Wang and Li [2023] have posted the optimization concerns on UL-based solvers. However, we call these claim into question by substantially outperforming heuristics DGA and RGA for the MIS on graphs with d = 20, 100, without training/historical instances D = {G µ } n µ=1 , as shown in Table 1 . See Appendix 6 for the results of solving all other Gsets, where consistently, CRA-PI-GNN provides better results as well.A comparison of the sampling-based solvers, RL-based solvers, SL-based solvers, Gurobi, and MIS-specific solvers is presented in Appendix F.2.We also compared the learning curve between PI-GNN and CRA-PI-GNN solver to confirmed that the CRA-PI-GNN solver does not become trapped in the plateau, p N = 0 N , as discussed in Section 3.1. Fig. 5 shows the dynamics of the cost functions for the MIS problems with N = 10,000 across d = 3, 5, 20, 100. Across all degrees, CRA-PI-GNN solver achieves a better solution with fewer epochs than PI-GNN solver. Specifically, PI-GNN solver becomes significantly slower due to getting trapped in the plateau even for graphs with low degrees, such as d = 3, 5. In contrast, CRA-PI-GNN solver can effectively escape from plateaus through the smoothing and automatic rounding of the penalty term when the negative parameter γ > 0. Computational scaling We next evaluate the computational scaling of the CRA-PI-GNN solver for MIS problems with large-scale RRGs with a node degree of 100 in Fig. 4 , following previous studies [Schuetz et al., 2022a, Wang and Li, 2023] . Fig. 4 demonstrated a moderate super-linear scaling of the total computational time, approximately ∼ N 1.4 for GCN and ∼ N 1.7 for GraphSage. This performance is nearly identical to that of the PI-GNN solver [Schuetz et al., 2022a] for problems on RRGs with lower degrees. It is important note that the runtimes of CRA-PI-GNN solver heavily depend on the optimizer for GNNs and annealing rate ε; thus this scaling remains largely unchanged for problems other than the MIS on 100 RRG. Additionally, CRA demonstrate that the runtime remains nearly constant as graph order and density increase, indicating effective scalability with denser graphs which is presented in Appendix F.2.Degree dependency of solutions We first compare the performances of PI-GNN and CRA-PI-GNN solvers with GCV, following Schuetz et al. [2022a] . Fig. 3 shows the cut ratio ν d as a function of the degree d compared to the theoretical upper bound [Parisi, 1980 , Dembo et al., 2017] . Across all degrees d, CRA-PI-GNN solver also outperforms PI-GNN solver, approaching the theoretical upper bound. In contrast, PI-GNN solver fails to find valid solutions for d > 20 as with the case of the MIS problems in Section 5.2. [Ye, 2003] , which is commonly used to evaluate MaxCut algorithms. Here, we provide benchmark results for seven distinct graphs with thousands of nodes, including Erdös-Renyi graphs with uniform edge probability, graphs in which the connectivity decays gradually from node 1 to N , 4-regular toroidal graphs, and a very large Gset instance with N = 10,000 nodes. Table 2 shows, across all problems, CRA-PI-GNN solver outperforms both the PI-GNN, RUN-CSP solvers and other greedy algorithm. See Appendix 6 for the results of solving all other Gsets, where consistently, CRA-PI-GNN provides better results as well. To evaluate the applicability of the CRA-PI-GNN solver to more practical problems not on graphs, we conducted experiments on DBM problems [Ferber et al., 2020 , Mulamba et al., 2020 , Mandi et al., 2022] ; refer to Appendix E.1 for details. This problems consists of 27 distinct instances with varying properties, and each instance comprises 100 nodes representing scientific publications, divided into two groups of 50 nodes N 1 and N 2 . The optimization is formulated as follows:l(x; C, M, λ) = -ij C ij x ij + λ 1 i ReLU j x ij -1 + λ 2 j ReLU i x ij -1 + λ 3 ReLU p ij x ij -ij M ij x ij + λ 4 ReLU q ij x ij -ij (1 -M ij )x ij ,where C ∈ R N1×N2 represents the likelihood of a link between each pair of nodes, an indicator M ij is set to 0 if article i and j share the same subject field (1 otherwise) ∀i ∈ N 1 , and j ∈ N 2 . The parameters p, q ∈ [0, 1] represent the probability of pairs sharing their field and of unrelated pairs, respectively. As in Mandi et al. [2022] , we explore two variations of this problem, with p = q = being 25% and 5%, respectively, and these variations are referred to as Matching-1 and Matching-2, respectively. In this experiment, we set λ 1 = λ 2 = 10 and λ 3 = λ 4 = 25. Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: This paper is accompanied by the code, along with detailed instructions to reproduce the experimental results. This will ensure that other researchers can faithfully reproduce the findings reported in the paper.Guidelines:• The answer NA means that paper does not include experiments requiring code. Justification: The research conforms to the NeurIPS Code of Ethics, ensuring responsible conduct and consideration of ethical implications.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 circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
conclusion: 
 This study proposes CRA strategy to address the both optimization and rounding issue in UL-based solvers. CRA strategy introduces a penalty term that dynamically shifts from prioritizing continuous solutions, where the non-convexity of the objective function is effectively smoothed, to enforcing discreteness, thereby eliminating artificial rounding. Experimental results demonstrate that CRA-PI-GNN solver significantly outperforms both PI-GNN solver and greedy algorithms across various complex CO problems, including MIS, MaxCut, and DBM problems. CRA approach not only enhances solution quality but also accelerates the learning process.Limitation In this numerical experiments, most hyperparameters were fixed to their default values, as outlined in Section 5.1, with minimal tuning. However, tuning may be necessary for certain problems or to further enhance performance.The theorem holds for the cross entropy penalty given byΦ(p) = N i=1 (p i log(p i ) + (1 -p i ) log(1 -p i )) (7)in the UL-based solver using data or history [Sun et al., 2022 , Sanokowski et al., 2024] EQUATIONC Generalization of CRA 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes]Justification: The paper discusses the potential positive impacts of the proposed CRA strategy, such as improving the efficiency of solving CO problems, as well as the potential negative impacts.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 or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness 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 a 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 or 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 a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper does not involve the release of high-risk data or models.Guidelines:• The answer NA means that the paper poses no 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 is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: 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?Answer: [Yes] Justification: The paper does not use any existing assets, so this question is not applicable.Guidelines:• The answer NA means that the paper does not 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 a 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. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • 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 are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: The paper does not introduce new assets.Guidelines:• The answer NA means that the paper does not release new 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. 
crowdsourcing 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 as details about compensation (if any)?Answer: [NA] Justification: The paper does not involve crowdsourcing or 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 contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to 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. 
 
 Appendix: This section generalize the penalty term Φ(θ; C) introduced for binary variables to K-Potts variables. K-Potts variable is the Kronecker delta δ(x i , x j ) which equas one whenever x i = x j and zero otherwise and a decision variable takes on K different values, ∀i ∈ [N ], x i = 1, . . . , K. For example, graph K-coloring problems can be expressed asEQUATIONFor Potts variables, the output of the GNN is h L θ = P (θ; C) ∈ R N ×K and the penalty term can be generalized as follows.EQUATIONWe describe the details of the GNN architectures used in all numerical experiments. The first convolutional layer takes H 0 -dimensional node embedding vectors, h 0 θ for each node, as input, yielding H 1 -dimensional feature vectors h 1 θ . Then, the ReLU function is applied as a componentwise nonlinear transformation. The second convolutional layer takes the H 1 -dimensional feature vectors, h 1 θ , as input, producing a H (2) -dimensional vector h 2 θ . Finally, a sigmoid function is applied to the H (2) -dimensional vector h 2 θ , and the output is the soft solution p θ ∈ [0, 1] N . As in [Schuetz et al., 2022a] , we set H 0 = int(N 0.8 ) or , H 1 = int(N 0.8 /2) and H 2 = 1 for both GCN and GraphSAGE .For all numerical experiments, we use the AdamW [Kingma and Ba, 2014] optimizer with a learning rate of η = 10 -4 and weight decay of 10 -2 , and we train the GNNs for up to 10 5 epochs with early stopping set to the absolute tolerance 10 -5 and patience 10 3 . As discussed in Schuetz et al. [2022a] , the GNNs are initialized with five different random seeds for a single instance because the results are dependent on the initial values of the trainable parameters; thus selecting the best solution.MIS problems There are some theoretical results for MIS problems on RRGs with the node degree set to d, where each node is connected to exactly d other nodes. The MIS problem is a fundamental NP-hard problem [Karp, 2010] defined as follows. Given an undirected graph G(V, E), an independent set (IS) is a subset of nodes I ∈ V where any two nodes in the set are not adjacent. The MIS problem attempts to find the largest IS, which is denoted I * . In this study, ρ denotes the IS density, where ρ = |I|/|V |. To formulate the problem, a binary variable x i is assigned to each node i ∈ V . Then the MIS problem is formulated as follows:EQUATIONwhere the first term attempts to maximize the number of nodes assigned 1, and the second term penalizes the adjacent nodes marked 1 according to the penalty parameter λ. In our numerical experiments, we set λ = 2, following Schuetz et al.[2022a], no violation is observed as in [Schuetz et al., 2022a] . First, for every d, a specific value ρ * d , which is dependent on only the degree d, exists such that the independent set density |I * |/|V | converges to ρ * d with a high probability as N approaches infinity [Bayati et al., 2010] . Second, a statistical mechanical analysis provides the typical MIS density ρ Theory d , as shown in Fig. 2 , and we clarify that for d > 16, the solution space of I undergoes a clustering transition, which is associated with hardness in sampling [Barbier et al., 2013] because the clustering is likely to create relevant barriers that affect any algorithm searching for the MIS I * . Finally, the hardness is supported by analytical results in a large d limit, which indicates that, while the maximum independent set density is known to have density ρ * d→∞ = 2 log(d)/d, to the best of our knowledge, there is no known algorithm that can find an independent set density exceeding ρ alg d→∞ = log(d)/d [Coja-Oghlan and Efthymiou, 2015].The MaxCut problem is also a fundamental NP-hard problem [Karp, 2010] with practical application in machine scheduling [Alidaee et al., 1994] , image recognition [Neven et al., 2008] and electronic circuit layout design [Deza and Laurent, 1994] . The MaxCut problem is also a fundamental NP-hard problem [Karp, 2010] Given an graph G = (V, E), a cut set C ∈ E is defined as a subset of the edge set between the node sets dividing(V 1 , V 2 | V 1 ∪ V 2 = V, V 1 ∩ V 2 = ∅).The MaxCut problems aim to find the maximum cut set, denoted C * . Here, the cut ratio is defined as ν = |C|/|V|, where |C| is the cardinality of the cut set. To formulate this problem, each node is assigned a binary variable, where x i = 1 indicates that node i belongs to V 1 , and x i = 0 indicates that the node belongs to V 2 . Here, x i + x j -2x i x j = 1 holds if the edge (i, j) ∈ C. As a result, we obtain the following:EQUATIONThis study has also focused on the MaxCut problems on d-RRGs, for which several theoretical results have been established. Specifically, for each d, the maximum cut ratio is given by ν* d ≈ d/4 + P * d/4 + O( √ d), where P * = 0.7632 . . . with a high probability as N approaches infinity [Parisi, 1980 , Dembo et al., 2017] . Thus, we take ν UB d = d/4 + P * d/4 as an upper bound for the maximum cut ratio in the large n limit.DBM problems Here, the topologies are taken from the Cora citation network [Sen et al., 2008] , where each node has 1,433 bag-of-words features, and each edge represents likelihood, as predicted by a machine learning model. Mandi et al. [2022] focused on disjoint topologies within the given topology, and they created 27 distinct instances with varying properties. Each instance comprises 100 nodes representing scientific publications, divided into two groups of 50 nodes N 1 and N 2 . The optimization task is to find the maximum matching, where diversity constraints ensure connections among papers in the same field and between papers of different fields. This is formulated using a penalty method as follows.EQUATION)where C ∈ R N1×N2 represents the likelihood of a link between each pair of nodes, an indicator M ij is set to 0 if article i and j share the same subject field (1 otherwise) ∀i ∈ N 1 , and j ∈ N 2 . The parameters p, q ∈ [0, 1] represent the probability of pairs sharing their field and of unrelated pairs, respectively. As in Mandi et al. [2022] , we explore two variations of this problem, with p = q = being 25% and 5%, respectively, and these variations are referred to as Matching-1 and Matching-2, respectively. In this experiment, we set λ 1 = λ 2 = 10 and λ 3 = λ 4 = 25.A GNN [Gilmer et al., 2017 , Scarselli et al., 2008] is a specialized neural network for representation learning of graph-structured data. GNNs learn a vectorial representation of each node in two steps, i.e., the aggregate and combine steps. The aggregate step employs a permutation-invariant function to generate an aggregated node feature, and in the combine step, the aggregated node feature is passed through a trainable layer to generate a node embedding, known as "message passing" or the "readout phase." Formally, for a given graph G = (V, E), where each node feature h 0 v ∈ R H0 is attached to each node v ∈ V , the GNN updates the following two steps iteratively. First, the aggregate step at each l-th layer is defined as follows:EQUATION)where the neighborhood of v ∈ V is denotedN v = {u ∈ V | (v, u) ∈ E}, h l-1 uis the node feature of the neighborhood, and a l v is the aggregated node feature of the neighborhood. Then, the combined step at each l-th layer is defined as follows:EQUATIONwhere h l v ∈ R H l denotes the node representation at the l-th layer. Here, the hyperparameters for the total number of layers L and the intermediate vector dimension N l are determined empirically. Although numerous implementations of GNN architectures have been proposed to date, the most basic and widely used architecture is the GCN [Scarselli et al., 2008] , which is given as follows:EQUATIONwhere W l and B l are trainable parameters, |N (v)| serves as a normalization factor, and σ : R H l → R H l is some component-wise nonlinear activation function, e.g., the sigmoid or ReLU function.F.1 Numerical validation of practical issues presented in Section 3.1In this sectioin, we will examine the issue (I) with continuous relaxations and the issue (II), the difficulties of optimization, as pointed out by previous studies [Wang and Li, 2023, Angelini and Ricci-Tersenghi, 2023] , in the NP-hard problems of MIS and the MaxCut problem. Therefere, we conducted numerical experiments using the PI-GNN solver for MIS and MaxCut problems on RRGs with higher degrees. To ensure the experimental impartiality, we adhered to the original settings of the PI-GNN solver [Schuetz et al., 2022b] . Refer to Section E for the detailed experimental settings. Fig. 7 (top) shows the solutions obtained by the PI-GNN solver as a function of the degree d for the MIS and MaxCut problems with varying system sizes N . These results indicate that finding independent and cut sets becomes unfeasible as the RRG becomes denser. In addition, to clarify the reasons for these failures, we analyzed the dynamics of the cost function for MIS problems with N = 10,000, with a specific focus on a graph with degrees d = 5 and d = 20, as depicted in Fig. 7 (bottom). For the d = 5 case, the cost function goes over the plateau of l(θ; G, λ) = 0 with p θ (G) = 0 N , as investigated in the histogram, eventually yielding a solution comparable to those presented by Schuetz et al. [2022a] . Conversely, in the d = 20 case, the cost function remains stagnant on the plateau of l(θ; G, λ) = 0 with p θ (G) = 0 N , thereby failing to find any independent nodes. Interpreting this phenomenon, we hypothesize that the representation capacity of the GNN is sufficiently large, leading us to consider the optimization of LMIS (θ; G, λ) and LMaxCut (θ; G) as a variational optimization problem relative to p θ . In this case, p * θ = 0 N satisfies the first-order variational optimality conditions δ lMIS /δp θ | p θ =p * = δ lMaxCut /δp θ | p θ =p * = 0 N , which implies a potential reason for absorption into the plateau. However, this does not reveal the conditions for the convergence to the fixed point p * during the early learning stage or the condition to escape from the fixed point p * . Thus, an extensive theoretical evaluation through stability analysis remains an important topic for future work.In summary, UL-based solver, minimizing θ can be challenging and unstable. In particular, the PI-GNN solver, which is one of the UL-based solvers employing GNNs, fails to optimize θ due to a local solution in complex CO problems on relatively dense graphs where the performance of greedy algorithms worsens. This issues can be potential bottleneck for more practical and relatively dense problems, making it challenging to employ the PI-GNN solver confidently.We evaluate our method using the MIS benchmark dataset from recent studies [Goshvadi et al., 2023 , Qiu et al., 2022] , which includes graphs from SATLIB [Hoos and Stützle, 2000] and Erdős-Rényi graphs (ERGs) of varying sizes. Following Sun et al. [2023] , our test set consists of 500 SATLIB graphs, each containing between 403 and 449 clauses with up to 1,347 nodes and 5,978 edges, 128 ERGs with 700 to 800 nodes each, and 16 ERGs with 9,000 to 11,000 nodes each. We conducted numerical experiments on PQQA using four different configurations: parallel runs with S = 100 or S = 1000 and shorter steps (3000 steps) or longer steps (30000 steps), similar to the approach in iSCO [Sun et al., 2023] . Table F .2 presents the solution quality and runtime results. The results show that CRA, which optimizes the relaxed variables as an optimization of GNN parameters, takes extra time for smaller ER-[700-800] instances due to the smaller number of decision variables. However, for larger instances, CRA achieves results comparable to iSCO. Although limited space makes it difficult to present other benchmark results employed by iSCO, such as MaxCut and MaxClique, numerical experiments on these benchmarks also show that CRA is less effective for small problems. However, for larger problems, the results are comparable to or slightly inferior to those of iSCO.We also investigated the relationship between the order of the graph and the solving time of the solver, and the results are shown in Table F .2 and F.2. The results demonstrate that the runtime remains nearly constant as graph order and density increase, indicating effective scalability with denser graphs.We conducted experiments across the additional GSET collection to further validate that including CRA enhances PI-GNN results beyond previously achievable in Table 6 .Table 3 : ApR and runtime are evaluated on three benchmarks provided by DIMES [Qiu et al., 2022] . The ApR is assessed relative to the results obtained by KaMIS. Runtime is reported as the total clock time, denoted in seconds (s), minutes (m), or hours (h). The runtime and solution quality are sourced from iSCO [Sun et al., 2023] . The baselines include solvers from the Operations Research (OR) community, as well as data-driven approaches utilizing Reinforcement Learning (RL), Supervised Learning (SL) combined with Tree Search (TS), Greedy decoding (G), or sampling (S). Methods that fail to produce results within 10 times the time limit of DIMES are marked as N/A. We conducted additional experiments on several TSP problems from the TSPLIB dataset1 , presenting results that illustrate the α-dependency.Experiments calculated the ApR as the ratio of the optimal value to the CRA result, with the ApR representing the average and standard deviation over 3 seeds. The "-" symbol in PI-GNN denotes cases where most variables are continuous values and where no solution satisfying the constraint was found within the maximum number of epochs. The same GNN and optimizer settings were used as in the main text experiments in Section 5.1. Table F .4 shows that the CRA approach yielded solutions with an ApR exceeding 0.9 across various instances. Notably, for the burma14 problem, our method identified the global optimal solution (3,323) multiple times. However, the optimal value can vary based on the specific GNN architecture and problem structure, indicating that a more comprehensive ablation study could provide valuable insights in future work.We conducted an ablation study focusing on the initial scheduling value γ(0) and scheduling rate ε. This numerical experiment was conducted under the configuration described in Section 5 and E. Fig. 8 shows the IS density of N = 10000 MIS problems on a 100-RRG as a function of the initial scheduling value γ(0) and the scheduling rate ε using the CRA-PI-GNN with both GraphSage and GCV. As can be seen, smaller initial scheduling γ(0) and scheduling rate ε values typically yield better solutions. However, the convergence time increases progressively as the initial scheduling γ(0) and scheduling rate ε values become smaller. In addition, GraphSage consistently produces better solutions even with relatively larger initial scheduling γ(0) and scheduling rate ε values, which implies that the GNN architecture influences both the solution quality and the effective regions of the initial scheduling γ(0) and scheduling rate ε values for the annealing process.Next, we investigated the effect of varying the curvature α in Eq. 5. Numerical experiments were performed on MIS problems with 10,000 nodes and the degrees of 5 and 20, as well as MaxCut problems with 10,000 nodes and the degrees of 5 and 35. The GraphSAGE architecture was employed, with other parameters set as Section 5 adn E. As shown in Fig. 9 , we observed that α = 2 consistently yielded the best scores across these problems.dataset, or 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., 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 all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to 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 be possible for other researchers to have some path to reproducing or verifying the results.http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/ 


Decision: While the original submission exhibited some minor flaws (such as writing issues, lack of comparison on a routing problem like the TSP, more discussion in the experiments section),  the authors have already dealt with one (TSP) and committed to fixing the others in a camera ready version. I expect them to follow through on this commitment. 

This paper presents a nice, simple idea, that works very well. I further emphasize that the writing quality is not good enough, but a quick check with a professional proofreader or LLM ought to fix that easily. Reviewer BuGW also provides a number of small corrections. The experiments, especially with the addition of the TSP in the review phase, provide a convincing argument that the simple annealing mechanism introduced is effective at solving a simple, but diverse set of combinatorial optimization problems. The unsupervised learning technique is a fascinating way forward that stands in contrast to many of the other learning-heavy proposals for neural combinatorial optimization.

Reviews: review: paper_summary: paper_summary: This paper presents a heuristic method for producing solutions to combinatorial optimization problems, which is based around solving a continuous relaxation of the problem. The main focus of the paper is on an additional penalty term to add to the objective of this relaxation which aims to reward solutions that are closer to satisfying the integrality constraints on the decision variables.
strength_weakness: strength_weakness: The paper is very dense and hard to follow, with little context provided to the reader. The method presented and evaluated in the computational study is ultimately an extension of the "PI-GNN" solver, but this fact is oddly kind of buried, with only an indirect reference in the introduction ("the solver that applies the CRA to the PI-GNN solver is referred to as the CRA-PI-GNN solver", with no indication that this is a main takeaway from the work), and then again at the end of Section 3.2. The paper does not explain in detail or formality what the PI-GNN solver is or how it works (how, specifically, does the CGA actually hook into PI-GNN?), and so a reader without prior familiarity cannot really understand or assess the new contributions laid out in Section 3. Ultimately, I do not feel confident that I can understand, and thus evaluate, the contributions proposed in the paper.
questions: questions: The main contribution of the paper is an additional penalization term to induce solutions that are feasible w.r.t. the binary constraints on the decision variables, but I do not see explicit discussion in the Experiments section about the feasibility of the solutions produced (wr.t. both the integrality constraints and the other equality/inequality constraints). Are all of the solutions used in the computational study feasible for the "true" problem? If there are numerical tolerances used to "fudge" exact feasibility, what are these tolerance values?
limitations: limitations: N/A
paper_summary: paper_summary: This paper aims to tackle shortcomings of the existing unsupervised learning-based solvers for combinatorial optimization, namely the local optima issue and the rounding issue. It proposes a novel technique called continuous relaxation annealing (CRA) strategy which introduces an additional penalty term to smooth the non-convexity of the objective function. This strategy is empirically shown to not only enhance the solution quality but also accelerate the learning process.
strength_weakness: strength_weakness: 1. My main concern lies in the technical contribution from this paper. The whole framework and empirical evaluation is built upon PI-GNN, which makes the observation and conclusion from this paper not generalizable.
2. I feel the research from this paper is kind of out-of-date. Check https://openreview.net/forum?id=ZMP0Bki9aK for SOTA results on the CO problems considered in this paper. In fact, [1] also mention that simulated annealing would perform better than GNN, but only greedy methods are used as baselines.


[1] Maria Chiara Angelini and Federico Ricci-Tersenghi. Modern graph neural networks do worse than classical greedy algorithms in solving combinatorial optimization problems like maximum independent set. Nature Machine Intelligence, 5(1):29–31, 2023.
questions: questions: N/A
limitations: limitations: Yes.
paper_summary: paper_summary: The proposed approach is an optimization method for each graph over GNN parameters where each output corresponds to the likelihood of the node belonging to the solution. The objective function consists of a penalty term along with a parameter scheduled to control the non-convexity of the objective.
strength_weakness: strength_weakness: [Major Comments]

1- The need to solve graph-based NP-hard problems that are originally formulated as ILPs stems from the unscalability of these solvers. For example, the scalability of the MIS problem depends on the number of nodes and the number of edges in the graph. This needs to be the motivation instead of the issues encountered in UL-based solvers.

2- While the proposed approach does not require training data (labeled or unlabeled), there are several hyper-parameters. Tuning these hyper parameters is a challenge. Further discussion is needed here. 

3- Getting trapped in local minima is not only the case in GNNs or PI-GNNs. It exists for any continuous relaxation of Problem 1. This is due to the non-convexity inherited in these formulations. For example, if we re-write Problem (3) in matrix form, we can see that the objective has a constant hessian equal to the adjacency matrix of the graph. If the magnitudes and signs of the eigen values vary significantly, then this indicates possible positive and negative curvatures in the loss landscape. Replacing x in Problem (3) with the output of a GNN does not guarantee changes. Although it may be possible that it will make some local minima avoidable by adaptive optimizers (such as ADAM), there is a possibility that this type of overparameterization would create unwanted local minima that do not result in any feasible solutions. Theoretically analyzing this is very complicated due to the use of a GNN. However, empirical investigation can be used to better motivate and understand the proposed approach. 

4- Similar to the previous point, rounding issues existed even before GNNs. See the SDP relaxations of MIS [1] and MaxCut [2] and how their dependence on rounding techniques (e.g. spectral clustering [3]) often fails to obtain optimal solutions. Rewriting is needed here. 

5- The Stationary point p* = 0_n was not discussed in Section 3.1. Furthermore, in line 208, it is 0_n, whereas in line 2016, it is 0_N. 

6- How was the GW approximation applied for the MaxCut problem? This approximation requires a normalized random vector drawn from the standard Gaussian distribution. How many samples were drawn? Given the SDP solution, one can simply draw multiple samples and pick the best where the only requirement is matrix-vector multiplication. This runs extremely fast with (i) no parameters of a NW, and (ii) no hyper-parameters to tune. The scenarios where such approaches fail need to be the motivation to propose the over-parameterized approach with convex annealing. 

7- Missing many "data-independent" baselines (methods that do not require pre-trained models (such as DIFUSCO [4]) or training data such as RL-based solvers (LwD [5])) for comparison such as ILP solvers (Gurobi, CPLEX, or CP-SAT [6]), sampling methods such as iSCO [7], SOTA heuristics such as ReduMIS [8], and differentiable solvers such as [9]. 

8- Why does the paper only consider d-regular graphs? How about the performance on other graphs? How does the run-time of this method scale in terms of the graph order and density? This is a major limitation of this work.

[Minor Comments]

1- What is script C in line 83?

2- What is I and J in the equation after line 86?

3- "nural" in line 106.

4- Cite equation 3. An example is [10]

5- Paragraph 149 to 151 is ill-sentenced.

6- “Indeed” in line 232.

7- Cite Potts variable optimization.

8- This study “employs” in line 237.

9- "are" in Appendix F.3 in line 247.

[References]

[1] On the shannon capacity of a graph. IEEE TIT, 1979.

[2] Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming. JACM, 1995.

[3] A tutorial on spectral clustering. Springer, 2007.

[4] Difusco: Graph-based diffusion solvers for combinatorial optimization. NeurIPS, 2023.

[5] Learning What to Defer for Maximum Independent Sets. ICML, 2020.

[6] https://developers.google.com/optimization

[7] Revisiting sampling for combinatorial optimization. ICML, 2023.

[8] A differentiable approach to the maximum independent set problem using dataless neural networks. Neural Networks, 2022.

[9] A branch and bound algorithm for the maximum clique problem. Computers & operations research, 1992.
questions: questions: See Weaknesses.
limitations: limitations: See Weaknesses.
paper_summary: paper_summary: This article finds that the existing UL-solvers will trap into local optima and face rounding issues. This study proposes a continuous relaxation annealing (CRA) strategy and an auxiliary function to facilitate training.
strength_weakness: strength_weakness: There are no major drawbacks in this article. There should be more reviews on neural combinatorial optimization solvers that apply annealing ideas as well (such as [1] provides annealing on the distance matrix of TSP).

[1] Lin, Xi, et al. ""Continuation path learning for homotopy optimization."" International Conference on Machine Learning. PMLR, 2023.
questions: questions: 1. Figure 8 provides parameter analysis for the N=10000 MIS problem. Can parameter sensitivity analysis be provided for other CO problems to demonstrate that CRA can be widely applied to general CO problems without special parameter design?
2. I am also concerned about the convergence speed under different parameter settings. Could you provide it as a function of the initial
scheduling and scheduling rate?
3. Can CRA be applied to routing problems such as TSP? If I understand correctly, the current $\phi$ function will have very small $p$-values in solving TSP, which may probably lead to a failure situation of CRA.
limitations: limitations: Well discussed.

Comments: title: Reply to qMFs
Reply to Reviewer J6G3
Reply to Reviewer qMFs
Response to Authors
Further response to reviewer BuGW
Thank you for the efforts. Response to Authors
comment: Thank you very much for your continued engagement and for addressing our responses. We greatly appreciate your acknowledgment of CRA's generalizability and the significance of our work.

> Firstly, I'm confused why authors feel there is any essential difference between "a first-order Taylor expansion of the objective function" used in iSCO and the "gradient of the objective function" used in CRA. 

We apologize for our unclear response. There might be a need for clarification regarding the iSCO [1]. We understand that **iSCO, which employs discrete Langevin Monte Carlo, is not a gradient descent algorithm based on continuous relaxation.** 
Instead, it generalizes Langevin dynamics (gradient flow) from continuous probability distributions to discrete ones. 
iSCO changes the temperature while numerically simulating the Equation from Theorem 3.1, characterized by Equation (23) on Sun 2022 [2].

Executing discrete Langevin dynamics without approximation requires calculating the difference between the objective functions $f(y)$ and $f(x)$ for the next transition candidates, as mentioned in Section 3.2.1 of Sun 2023 [1]. This approach is noted to be computationally intensive, as Sun et al. themselves acknowledge in Section 3.2.1. 
**The need to compute the difference between objective values rather than gradients indicates that iSCO is not based on continuous relaxation**. Although they approximate this difference with a Taylor expansion to expedite the simulation, it remains unclear how effective this approximation is for more complex problems and how to related to our method.

> Then it leads to my second question. Clearly, iSCO owns more parameters than CRA when the graph size is even larger, then why does CRA outperform iSCO on large instances if the over-parameterization is its advantage against iSCO? My personal answer here is that iSCO is optimized on a larger parameter space and it converges slowly on large instances.

As explained above, **iSCO uses the gradient of the objective function to approximate and efficiently simulate original discrete Langevin dynamics without using continuous relaxation**. 
On the other hand, since CRA uses continuous relaxation, it is easy to implement gradient descent, and optimizers with good convergence properties, such as AdamW, can be used. Indeed, we used AdamW in our numerical experiments. Also, as reviewer BuGW pointed out, over-parameterizing the decision variables using GNNs may smooth the loss landscape, making it more suitable for optimization and speeding up the convergence speed. 
However, it is difficult to investigate this hypothesis at present theoretically, and we consider that a more detailed investigation is future work. Additionally, it is crucial to investigate iSCO's implementation details, including its specific code-level parallelization methods. Unfortunately, as the iSCO code has yet to be shared, it has been challenging to delve into these details during the discussion period.  In the revised version, we will add these differences between iSCO and CRA to the related work.

> I note that all the theoretical analyses are on the output parameters $p_{i}$, isn't this type of analyses more related to gradient-based simulated annealing methods whose parameters are exactly $p_{i}$ here?

> I do question the significance of this work on the Type II solvers. 

Our analysis is a theoretical analysis of the output parameters $p_{i}$. 
As you point out, there is an indirect relationship with gradient-based simulated annealing for $p_{i}$, but gradient-based simulated annealing is a method that performs Langevin Monte Carlo while lowering the temperature, and our method is a gradient descent method that controls the term $\Phi(p)$ that controls the continuity and discreteness.
**Also, the method of simply performing continuous relaxation and gradient-based simulated annealing for $p_{i}$ cannot guarantee the original CO problems because they have the rounding issue, as pointed out in our main text.** 

We observed superior results compared to DIFUSCO for larger instances. As a heuristic approach, obtaining good solutions for problems too large for an exact solution method like Gurobi is vital. Also, as Reviewer BuGW pointed out, a more comprehensive comparison is needed between CRA and iSCO, including the setting of hyperparameters. We think that investigating which method is superior in a more quantitative way is future work (While we intended to conduct a detailed comparison with iSCO on larger problems such as d-regular random graphs during the rebuttal period, the unavailability of the iSCO code prevented us from completing this within the time limit).

Considering these points, we would sincerely appreciate it if you could provide a final evaluation.

- [1] Sun et al. Revisiting Sampling for Combinatorial Optimization. ICML 2023.
- [2] Sun et al. Discrete langevin samplers via wasserstein gradient flow. AISTATS 2023.
Thank you for your additional response. 

Generally, I'm satisfied with the technical part of this work, but I'm largely confused by the takeaways from this paper as it is different from my intuition. There could be bias from my personal research, so I promise to authors that I will continue the discussion with other reviewers to see if I take it incorrectly.

- Firstly, I'm confused why authors feel there is any essential difference between "a first-order Taylor expansion of the objective function" used in iSCO and the "gradient of the objective function" used in CRA. From my point of view, gradient is exactly an outcome of the first-order Taylor expansion. In iSCO, the parameter space spans the scores $p_i$ on all nodes/edges; and in CRA, the parameter space spans all the GNN parameters. 

- Then it leads to my second question. Clearly, iSCO owns more parameters than CRA when the graph size is even larger, then why does CRA outperform iSCO on large instances if the over-parameterization is its advantage against iSCO? My personal answer here is that iSCO is optimized on a larger parameter space and it converges slowly on large instances. 

The reason I stated that this research is out-of-date is because the method introduced by Wang (2023) is clearly worse than iSCO in both its performance and the interpretability. Given the recent advances in sampling-based method, either unsupervised (iSCO) or supervised (DFUSCO), I do question the significance of this work on the Type II solvers. I note that all the theoretical analyses are on the output parameters $p_i$, isn't this type of analyses more  related to gradient-based simulated annealing methods whose parameters are exactly $p_i$ here?

I do agree with authors that "CRA can be easily generalized to UL-based solvers (Type I) and other relaxation-based solvers", which seems to be a more significant topic to me. But currently, I feel this work is quite restricted to a narrow problem where other methods have already stood out.
Thank you for recognizing the contribution of our paper, and we appreciate your decision to raise the score.

Based on your insightful comments, we conducted additional experiments on several TSP problems from the TSPLIB dataset (http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/). Below, we present the results demonstrating the $p$-dependency of the solutions.

| Instance   | ApR ($p=2$)          | ApR ($p=4$)          | ApR ($p=6$)          | ApR ($p=8$)          | ApR (PI)            | Optimum |
|------------|----------------------|----------------------|----------------------|----------------------|---------------------|---------|
| burma14    | 0.91 ± 0.08          | 0.98 ± 0.10          | 0.97 ± 0.14          | 0.99 ± 0.06          | 0.736 ± 1.21        | 3,323   |
| ulysses22  | 0.89 ± 0.03          | 0.92 ± 0.02          | 0.88 ± 0.07          | 0.89 ± 0.05          | --                  | 7,013   |
| st70       | 0.96 ± 0.01          | 0.85 ± 0.03          | 0.88 ± 0.04          | 0.80 ± 0.02          | --                  | 675     |
| gr96       | 0.81 ± 0.05          | 0.82 ± 0.03          | 0.90 ± 0.05          | 0.86 ± 0.05          | --                  | 5,5209  |


In these experiments, we calculated the ApR as the ratio of the optimal value to the CRA result, with the ApR representing the average and standard deviation over 3 seeds. 
The "--" in PI-GNN indicates that most variables are continuous values and that the maximum epoch did not obtain a solution satisfying the constraint. 
We used the same GNN and optimizer as in our experiments of the main text.

The table shows that the CRA approach found solutions with an ApR exceeding 0.9 across various instances. Notably, for the "burma14" problem, our method found the global optimal solution, 3,323, several times. However, the optimal $p$ value may vary depending on the specific GNN architecture and problem structure, which suggests that a more comprehensive ablation study would be beneficial in future work.

We will include these results and additional problem instances in the appendix or main text of the revised version.
We hope that these enhancements, which demonstrate the broader applicability and robustness of our approach, will further highlight the contribution of our work. We would greatly appreciate it if these improvements could be taken into consideration during your final evaluation.
Thank you for your thoughtful feedback and continued engagement with our work. 

We apologize for any confusion caused by our previous explanation in the global (unified) response. 
The iSCO results presented in the *"Numerical Experiments"* of our global response are from the paper by Sun et al. (ICML 2023) [1]. 
We apologize for any inconvenience and kindly request that you review these results in the global response, where we demonstrated CRA's performance that is comparable to or slightly below the SOTA iSCO method.
As Reviewer BuGW noted, **"I acknowledge the additional experiments for MIS. The method clearly outperforms DIFUSCO while requiring no training data. However, the proposed approach still underperforms when compared to iSCO, a training-data-free sampling method where no GNNs are needed. However, I think that iSCO, while reporting significant results, is sensitive to the choice of the sampler and hyper-parameters. So overall, the CRA results in this table are fair."** This recognition of our comparison with iSCO supports the fairness of our results.

It is also important to highlight that iSCO [1] uses a first-order Taylor expansion of the objective function to approximate the essential $\Delta(x)$ term in the discrete Langevin dynamics, enabling efficient parallel processing on GPUs. 
However, this method may experience performance degradation when the first-order approximation is inadequate. The iSCO paper primarily tests problems where the approximation is valid or nearly valid, which should be considered, particularly in black-box optimization with complex surrogate models.

In contrast, our method does not rely on such approximations and can be applied as long as the gradient of the objective function, $\hat{l}(p; C, \lambda)$, is accessible. Indeed, Wang et al. (2022) [2] demonstrated the superiority of UL-based solvers in black-box optimization using complex surrogates with GNNs, leveraging the flexibility of UL-based solvers. 
Additionally, as Reviewer BuGW pointed out, **"I believe that over-parameterization and benefiting from the structure of a GNN are what makes the contribution valuable. To me, the benefit here is similar, in spirit, to the topic of lifting in optimization."** 
This suggests that solution performance could further improve depending on the GNN's over-parameterization and architecture.
We believe that the statement in the *"Main Contribution and Novelty"* section of the Global Response, "Our primary contribution is the introduction of CRA to overcome the limitations of UL-based solvers (Type II), as highlighted by Wang et al. (ICLR 2023)," is crucial to realizing these future possibilities.

To further strengthen our paper, we plan to include more detailed comparisons with iSCO in the revised version and additional discussions and results. We respectfully request that you reconsider your score based on these efforts.

- [1] Sun et al. Revisiting Sampling for Combinatorial Optimization. ICML 2023.
- [2] Wang et al. Unsupervised learning for combinatorial optimization with principled objective relaxation. NeurIPS 2022.
I sincerely appreciate authors' for your response. However, I believe the link included in my second weakness part refers to another paper [1] that
- achieves the SOTA performance
- "performs gradient-based optimization and thus efficiently utilizes GPUs and state-of-the-art optimizers" exactly as you claimed for your method.
- and also requires no supervision

I do not see a clear discussion here.

[1] Sun et al. Revisiting Sampling for Combinatorial Optimization. ICML 2023.
I would like to thank the authors for their response to my question and addressing most of my concerns. The results of the attached PDF make more sense now. Initially, I was looking at the increasing $d$ and $p$ for density as these parameters indicate the density of the these graphs. I am happy to further increase my score. 

Final remarks:

- While the run-time comparison with KaMIS in Table 2 is encouraging, I wouldn't call CRA solution quality comparable. KaMIS is a CPU-based C++ (or C) implementation. I wonder whether it can be designed to be faster with GPUs (just sharing a thought, not asking for more experiments). However, this method depends on MIS-specific graph reductions that would still be slower on denser graphs. 

- **Overall**, I feel that the authors should emphasize that their method is dataless (no training data is required) which, I believe, is a great advantage and should have been reflected in the title. I believe that over-parameterization and benefiting from the structure of a GNN are what makes the contribution valuable. To me, the benefit here is similar, in spirit, to the topic of lifting in optimization.
Thank you for raising your score, thoroughly reviewing it, and providing insightful comments. 
We also thank you for your insightful understanding of over-parameterization. 

> However, first, I find that the Apr in these results are low, and KaMIS should have been used to examine the solution quality. 

We apologize for any confusion caused by the results in our global response's PDF. 
**The results are based on IS density $\rho$ as defined in Line 99, not ApR.**
To clarify, we have included below a revised table where ApR is calculated, comparing our method against theoretical results [Barbier et al., 2023], as in Lines 299-302 in our experimental section.

| Problem               | $\mathrm{ApR}$ (CRA) | $\mathrm{ApR}$ (PI) | Time (CRA) | Time (PI)  |
|-----------------------|--------------|-------------|------------|------------|
| $\mathrm{RRG}(1{,}000, 10)$ | 0.95         | 0.78        | 108 (s)    | 98 (s)     |
| $\mathrm{RRG}(1{,}000, 20)$ | 0.95         | 0.56        | 103 (s)    | 92 (s)     |
| $\mathrm{RRG}(1{,}000, 30)$ | 0.94         | 0.00        | 102 (s)    | 88 (s)     |
| $\mathrm{RRG}(1{,}000, 40)$ | 0.93         | 0.00        | 101 (s)    | 82 (s)     |
| $\mathrm{RRG}(1{,}000, 50)$ | 0.92         | 0.00        | 102 (s)    | 82 (s)     |
| $\mathrm{RRG}(1{,}000, 60)$ | 0.91         | 0.00        | 101 (s)    | 91 (s)     |
| $\mathrm{RRG}(1{,}000, 70)$ | 0.91         | 0.00        | 101 (s)    | 86 (s)     |
| $\mathrm{RRG}(1{,}000, 80)$ | 0.91         | 0.00        | 102 (s)    | 93 (s)     |
| $\mathrm{RRG}(5{,}000, 10)$ | 0.93         | 0.77        | 436 (s)    | 287 (s)    |
| $\mathrm{RRG}(5{,}000, 20)$ | 0.95         | 0.74        | 413 (s)    | 280 (s)    |
| $\mathrm{RRG}(5{,}000, 30)$ | 0.95         | 0.00        | 419 (s)    | 283 (s)    |
| $\mathrm{RRG}(5{,}000, 40)$ | 0.94         | 0.00        | 429 (s)    | 293 (s)    |
| $\mathrm{RRG}(5{,}000, 50)$ | 0.94         | 0.00        | 418 (s)    | 324 (s)    |
| $\mathrm{RRG}(5{,}000, 60)$ | 0.93         | 0.00        | 321 (s)    | 302 (s)    |
| $\mathrm{RRG}(5{,}000, 70)$ | 0.92         | 0.00        | 321 (s)    | 325 (s)    |
| $\mathrm{RRG}(5{,}000, 80)$ | 0.92         | 0.000       | 330 (s)    | 305 (s)    |

As shown by these results, the ApR exceeds 0.9 for all values of $d$.

Additionally, we have conducted further comparisons with KaMIS for the Erdos–Renyi graph, focusing on runtime and ApR, which is evaluated by comparing our method against KaMIS. Due to time limitations, we constrained the running time for KaMIS, and the results below show the average ApRs and runtimes across five different random seeds.

| Problem                  | CRA($\mathrm{ApR}$) | PI($\mathrm{ApR}$) | Time (CRA) | Time (PI) | Time (KaMIS) |
|--------------------------|---------------------|--------------------|------------|-----------|--------------|
| $\mathrm{ERG}(1{,}000, 0.05)$ | 0.97                | 0.01               | 103 (s)    | 98 (s)    | 100 (s)      |
| $\mathrm{ERG}(1{,}000, 0.10)$ | 0.95                | 0.00               | 100 (s)    | 98 (s)    | 210 (s)      |
| $\mathrm{ERG}(1{,}000, 0.15)$ | 0.94                | 0.00               | 100 (s)    | 92 (s)    | 315 (s)      |
| $\mathrm{ERG}(1{,}000, 0.20)$ | 0.91                | 0.00               | 99 (s)     | 88 (s)    | 557 (s)      |
| $\mathrm{ERG}(1{,}000, 0.25)$ | 0.93                | 0.00               | 98 (s)     | 82 (s)    | 733 (s)      |
| $\mathrm{ERG}(1{,}000, 0.30)$ | 0.90                | 0.00               | 98 (s)     | 82 (s)    | 1000 (s)     |
| $\mathrm{ERG}(1{,}000, 0.35)$ | 0.92                | 0.00               | 99 (s)     | 91 (s)    | 1000 (s)     |
| $\mathrm{ERG}(1{,}000, 0.40)$ | 0.91                | 0.00               | 97 (s)     | 86 (s)    | 1000 (s)     |

These results demonstrate that our method performs comparably to KaMIS. The revised manuscript will include a more thorough comparison of larger node cases in the main text or appendices.

Given this evidence and the detailed comparisons, we believe we have addressed this significant concern. We respectfully request that you reconsider your score based on these results. We are committed to further enhancing our paper as suggested and will include a more detailed discussion and additional results in the revised version.
### Response to Authors: 

- The experiments are evaluated on RRG and ER. The run-time of the proposed method scales only with the size of the graph, which I find impressive. Your response "runtime remains nearly constant as graph order and density increase" is not accurate as it is clear that the run-time is ~4X more when $n$ increased. In subsequent versions of the paper, I highly recommend including KaMIS and ILP solvers in these tables, as they are both known to scale with density, not just graph size. However, first, I find that the Apr in these results are low, and KaMIS should have been used to examine the solution quality. Second, in general, I still think that the proposed method is under-evaluated. The authors should consider using several random graph generators from NetWorkX, especially since it is a dataless method. 

- Discussing the intuition of over-parameterization to smoothen the landscape is acknowledged. I highly recommend that the authors include their response to this comment in the main body of the paper.

- I acknowledge the additional experiments for MIS. The method clearly outperforms DIFUSCO while requiring no training data. However, the proposed approach still underperforms when compared to iSCO, a training-data-free sampling method where no GNNs are needed. However, I think that iSCO, while reporting significant results, is sensitive to the choice of the sampler and hyper-parameters. So overall, the CRA results in this table are fair.

- **Question**: Why are the CRA MIS size results in the "Numerical Experiments" global response competitive, but the results in the attached PDF are not (way below 1)?

- Larger graphs sizes (nodes) or harder instances can be any graph where ILPs and heuristics struggle in terms of either run-time or solution quality. For MIS, the authors can try mid-density (with respect to the complete graph) GNM graphs or ER (with $p > 0.4$) graphs. 

- Overall, I appreciate the dataless training approach and the theoretical support of the proposed method, and I encourage the authors to further investigate this direction. However, I feel that the paper could benefit from (i) a major revision in terms of writing (e.g., polishing the arguments about the non-convexity and rounding issues of general continuous relaxations of COPs) and (ii) a more comprehensive evaluation and explanation of the results. Therefore, I will only increase my score to 4.
Thank you for your response.
﻿
There are no major drawbacks to this article. I will increase the rating to 6 to support your acceptance. 
﻿
I still believe that exploring applications in routing problems such as TSP is important, and I hope to see this part in the rebuttal period and the following versions of your manuscript.