Controlling Continuous Relaxation for
Combinatorial Optimization

One approach, proposed by Karalias and Loukas (2020), aims to automatically learn effective heuristics from historical dataset instances ${\mathcal{D}}=\{C_{\mu}\}_{\mu=1}^{P}$ and then apply these learned heuristics to a new instance $C^{\ast}$, 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 ${\mathcal{D}}=(C_{\mu})$, sampled independently and identically from a distribution $P(C)$, the goal is to minimize the average loss function $\min_{{\bm{\theta}}}\sum_{\mu=1}^{P}l({\bm{\theta}};C_{\mu},{\bm{\lambda}})$. However, this method does not guarantee quality for a test instance, $C^{\ast}$. Even if the training instances ${\mathcal{D}}$ are extensive and the test instance $C$ follows $P(C)$, low average performance ${\mathbb{E}}_{C\sim P(C)}[\hat{l}(\theta;C)]$ may not guarantee a low $l(\theta;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.

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_{\theta}(G)$. Here, an $L$-layered GNN is trained to directly minimize $\hat{l}({\bm{\theta}};C,{\bm{\lambda}})$, taking as input a graph $G$ and the embedding vectors on its nodes, and outputting the relaxed solution ${\bm{p}}_{{\bm{\theta}}}(G)\in[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 ${\mathcal{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).

These requirements encompass a nonlinear cost function and interactions involving many-body interactions, extending beyond simple two-body interactions.

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 $\nicefrac{{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 $\nicefrac{{1}}{{2}}$, underscoring the limitations of current rounding techniques for UL-based solvers.

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, ${\mathcal{D}}=\{C_{\mu}\}_{1\leq\mu\leq 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 ${\mathcal{D}}$. However, in a practical setting, systematic methods for generating or collecting training datasets ${\mathcal{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, ${\mathcal{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, ${\bm{p}}_{{\bm{\theta}}}={\bm{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 ${\bm{p}}_{{\bm{\theta}}}={\bm{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.

To address these issues, we propose a penalty term to control the balance between discreteness and continuity in the relaxed variables, formulated as follows:

where $\gamma\in{\mathbb{R}}$ is a penalty parameter, and the even number $\alpha$ denote a curve rate. When $\gamma$ is negative, i.e., $\gamma<0$, the relaxed variables tend to favor the continuous space, smoothing the non-convexity of the objective function $\hat{l}({\bm{p}};C,{\bm{\lambda}})$ due to the convexity of the penalty term $\Phi({\bm{p}})$. In contrast, when $\gamma$ is positive, i.e., $\gamma>0$, the relaxed variables tend to favor discrete space, smoothing out the continuous solution into discrete solution. Formally, the following theorem holds as $\lambda$ approaches $\pm\infty$.

For the detailed proof, refer to Appendix B.1. Theorem 3.1 can be generalized for any convex function $\Phi({\bm{p}};C)$ that has a unique maximum at $\nicefrac{{{\bm{1}}_{N}}}{{2}}$ and achieves a global minimum for all ${\bm{p}}\in\{0,1\}^{N}$; an example is binary cross entropy $\Phi_{\mathrm{CE}}({\bm{p}})=\sum_{i=1}^{N}(p_{i}\log p_{i}+(1-p_{i})\log(1-p_{i}))$, introduced by Sun et al. (2022); Sanokowski et al. (2024) for the UL-based solvers (Type I). Additionally, the penalty term eliminates the stationary point ${\bm{p}}^{\ast}={\bm{0}}_{N}$ described in Section 3.1, preventing convergence to a plateau. For UL-based solvers, the penalty term is expressed as follows:

where $\Phi({\bm{\theta}};C)\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}\Phi({\bm{p}}_{\theta}(C))$. According to Theorem 3.1, setting a sufficiently large $\gamma$ value cases the relaxed variables to approach nearly discrete values. We can also generalize this penalty term $\Phi({\bm{\theta}};C)$, to Potts variables optimization, including coloring problems (Schuetz et al., 2022b), and mixed-integer optimization; refer to Appendix C.1.

We propose an annealing strategy that gradually anneals the penalty parameter $\gamma$ in Eq. (6). Initially, a negative gamma value, i.e., $\gamma<0$, is chosen to leverage the properties, facilitating broad exploration by smoothing the non-convexity of $\hat{l}({\bm{\theta}};C,{\bm{\lambda}})$ and eliminating the stationary point ${\bm{p}}^{\ast}={\bm{0}}_{N}$ to avoid the plateau, as discussed in Section 3.1. Subsequently, the penalty parameter $\gamma$ is gradually increased to a positive value, $\gamma>0$, with each update of the trainable parameters (one epoch), until the penalty term approaches zero, i.e., $\Phi({\bm{\theta}},C)\approx 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 $\Phi_{\mathrm{CE}}({\bm{p}})$ is infeasible for UL-based solvers when $\gamma>0$, as the gradient $\nicefrac{{\partial\Phi_{\mathrm{CE}}({\bm{p}})}}{{\partial p_{i}}}$ diverges to $\pm\infty$ at $0$ or $1$. In deed, when $\gamma=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 $\Phi_{\mathrm{CE}}({\bm{p}})$ makes the learning infeasible without additional techniques, such as gradient clipping. In contrast, the gradient of the penalty term in Eq. 5, $\nicefrac{{\partial\Phi({\bm{p}})}}{{\partial p_{i}}}$, is bounded within $[-2\alpha,2\alpha]$ for any $\gamma$, preventing the gradient divergence issue seen in $\Phi_{\mathrm{CE}}({\bm{p}})$. Additionally, by increasing $\alpha$, the absolute value of the gradient near $\nicefrac{{1}}{{2}}$ becomes smaller, allowing for control over the smoothing strength toward a discrete solution near $\nicefrac{{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., $\Phi({\bm{\theta}};C)\approx 0$. Various annealing schedules can be considered; in this study, we employ the following scheduling: $\gamma(\tau+1)\leftarrow\gamma(\tau)+\varepsilon$, where the scheduling rate $\varepsilon\in{\mathbb{R}}_{+}$ is a small constant, and $\tau$ 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 $\gamma(0)$ and the scheduling rate $\varepsilon$. Numerical experiments suggest that better solutions are obtained when $\gamma(0)$ is set to a small negative value and $\varepsilon$ is kept low. The ablation study are presented in Appendix F.5.

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 ${\mathcal{D}}=\{C_{\mu}\}_{\mu=1}^{P}$. 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.

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 $\eta=10^{-4}$ and weight decay of $10^{-2}$. The GNNs are trained for up to $5\times 10^{4}$ epochs with early stopping, which monitors the summarized loss function $\sum_{s=1}^{S}\hat{l}(P_{:,s})$ and the entropy term $\Phi(P;\gamma,\alpha)$ with tolerance $10^{-5}$ and patience $10^{3}$ epochs; Further details are provided in Appendix D.2. We set the initial scheduling value to $\gamma(0)=-20$ for the MIS and matching problems, and we set $\gamma(0)=-6$ for the MaxCut problems with the scheduling rate $\varepsilon=10^{-3}$ and curve rate $\alpha=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\in[N]$, we map $p_{\theta,i}$ into $0$ if $p_{\theta,i}\leq 0.5$ and $p_{\theta,i}$ into $1$ if $p_{\theta,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.

Following Karalias and Loukas (2020), We use the approximation ratio (ApR), formulated as $\mathrm{ApR}=f({\bm{x}};C)/f({\bm{x}}^{\ast};C)$, where ${\bm{x}}^{\ast}$ is optimal solution. For the MIS problems, we evaluate the ApRs using the theoretical optimal cost (Barbier et al., 2013) and the independent set density $\rho$ relative to the theoretical results. For the MaxCut problems on RRGs, we adopt the cut ratio $\nu$ 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.

First, we compare the performance of the PI-GNN and CRA-PI-GNN solvers using GCV, as in Schuetz et al. (2022a). Fig. 3 shows the independent set density $\rho_{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\geq 15$, as pointed by the previous studies (Angelini and Ricci-Tersenghi, 2023; 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 ${\mathcal{D}}=\{G^{\mu}\}_{\mu=1}^{n}$, 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, ${\bm{p}}_{N}={\bm{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 $\gamma>0$.

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. 1, following previous studies (Schuetz et al., 2022a; Wang and Li, 2023). Fig. 1 demonstrated a moderate super-linear scaling of the total computational time, approximately $\sim N^{1.4}$ for GCN and $\sim 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 $\varepsilon$; 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.

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 $\nu_{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.

Following Schuetz et al. (2022a), we next conducted additional experiments on standard MaxCut benchmark instances based on the publicly available Gset dataset (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.