Variational Annealing on Graphs for Combinatorial Optimization
Authors: Sebastian Sanokowski
Wilhelm Franz Berghammer
Sepp Hochreiter
Sebastian Lehner
Venue: NeurIPS 2023 poster (2023)
Keywords: Combinatorial Optimization
Entropy Regularization
Graph Neural Networks
Statistical Mechanics

Abstract: Several recent unsupervised learning methods use probabilistic approaches to solve combinatorial optimization (CO) problems based on the assumption of statistically independent solution variables. We demonstrate that this assumption imposes performance limitations in particular on difficult problem instances. Our results corroborate that an autoregressive approach which captures statistical dependencies among solution variables yields superior performance on many popular CO problems. We introduce Subgraph Tokenization in which the configuration of a set of solution variables is represented by a single token. This tokenization technique alleviates the drawback of the long sequential sampling procedure which is inherent to autoregressive methods without sacrificing expressivity. Importantly, we theoretically motivate an annealed entropy regularization and show empirically that it is essential for efficient and stable learning.

Variational Annealing on Graphs for Combinatorial Optimization
Sebastian Sanokowski
Wilhelm Berghammer
Sepp Hochreiter 1.2
Sebastian Lehner
ELLIS Unit Linz and LIT AI Lab
Institute for Machine Learning, Johannes Kepler University Linz, Austria sanokowski@m. jku.
Abstract
Several recent unsupervised learning methods use probabilistic approaches to solve combinatorial optimization (CO) problems based on the assumption statistically independent solution variables. We demonstrate that this assumption imposes performance limitations particular on difficult problem instances. Our results corroborate that an autoregressive approach which captures statistical dependen- cies among solution variables yields superior performance on many popular CO problems. We introduce subgraph tokenization in which the configuration of a set of solution variables is represented by a single token. This tokenization tech- nique alleviates the drawback of the long sequential sampling procedure which is inherent to autoregressive methods without sacrificing expressivity. Importantly, theoretically motivate annealed entropy regularization and show empirically that is essential for efficient and stable learning.
Introduction
Combinatorial optimization (CO) problems are central interest to a wide range of fields, including operations research, physics, and computational complexity (Papadimitriou and Steiglitz 1998). Since CO problems are typically NP-hard one would not expect that one can find the solutions to all arbitrarily large CO problem instances in polynomial time. However, these restrictions concerned with worst case scenarios over entire problem families. For instance finding the Minimum Independent Set of any graph. Consequently, it is not surprising that a large body of work is focused on the design of solution strategies that yield a particularly good performance on a restricted subset of all possible instances. There is growing interest in exploring deep learning approaches this restricted CO setting [Bello et al., (017). These methods aim to learn how to efficiently generate high quality solutions for CO problems and rely primarily learned algorithmic components rather than on -obllem-ppecffic, hand-crafted ones (Bengio et al. [2021]. Importantly, practicable methods should not rely on supervised training with solutions from other solvers since this limits the achievable performance of the learned model to that of the solver used to generate the training data (Yehuda) al. 2020 As a result the development of unsupervised methods [Karalias and Loukas 2020
Wang et al. 2022 Qiu et al. 2022 Karalias et al. [022 Min et al. 2022 Wang and Li| 2023 and of
Reinforcement Learning (RL) methods [Bello et al. 201 Khalil al. 2017 Kool al. 2019 Ahn
al.) 2020 Böther et al. 2022 Mazyavkina et al. 2021 for CO is active field of research. Tn this
work we are interested training models that generate solutions at inference contrast to approaches
OOur code is available at: https:://gtthub.comm/ http:///ggttub..ccmmmlljjjuuu/e CK-C/
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

that require problem instance specific training. A popular approach is to take representation of the CO problem instance as input to a deep learning model and to output the parameters of distribution over solutions that has its probability mass concentrated on regions with high solution quality. This probabilistic optimization approach to CO is used in the recent works of e.g- (Karalias and Loukas] 2020 Min et al. 2022 Qiu et al. 2022, Sun et al. 2022]. In these works the distribution over solutions is built on the assumption of mutual statistical independence of the individual solution parameters and can, consequently, be represented by product of Bernoulli distributions for the the individual solution variables. This simplification is typically refereed to as a mean-field approximation (MFA) and is
frequently used in various fields including e.g. statistical mechanics [Parisi 1988], Bayesian statistics
(Wainwright and Jordan 2008], and the analysis of neural networks [Mei et al. 2019]. However, the simplifying assumption of the MFA restricts the expressivity of the corresponding distributions which limits its applicability when the target distribution represents strong correlations [Jaakkola and Jordan [998]. However, replacing MFA with more expressive approaches is computationally expensive and requires careful regularization to ensure efficient and stable training. Based on these consideration our contributions can be summarized follows.
We demonstrate that the frequently used MFA imposes limits on the attainable solution quality for CO problems. Importantly, we introduce Variational Annealing on Graphs for Combinatorial Optimization ( (VAG-CO) a method that achieves new state-of-the of-the--rtt performances on popular CO problems by combining expressive autoregressive models with annealed entropy regularization. We provide a theoretical motivation for this entropy regularization via considerations of the sample complexity of related density estimation problems. By introducing sub-graph tokenization VAG-CO drastically reduces the number of necessary steps to generate solutions and therein significantly improves the training and inference efficiency.
2 Problem Description
Ising Formulation of CO. To introduce VAG-CO it is convenient to adopt a point of view on CO that is motivated by statistical mechanics. As shown in Lucas [014] many frequently encountered NP-complete CO problems can be reformulated in terms of a particular class of models known as Ising models which therefore allow the representation of different CO problem types in unified manner. In this context the CO problem is equivalent to finding the minimum of an energy function R. This function assigns an energy E(or) an N-tuple of discrete solution variables E: 1-1,1)) NN) E 1-1, which are often referred to as spins. The family of aforementioned a = (ai
Ising models is characterized by the following form of the energy function:
N F(a)= Σ Jijoro + Σ B
where the first term represents the interaction between spins o, E {-1,1) through couplings Juy E R. while the second term determines the magnitude B; E R of the contribution of an individual spin to the energy of the system. For brevity we will denote the parameters of the energy function (Ji, B,) simply E. Hence, E will represent CO problem instance while E(or) will denote the energy of a solution or under an energy function with parameters E. The Ising formulation for the four CO problem types studied in this work is given in Tab. [
Variational Learning. Next, we specify how we approach the learning problem of approximating the optimal solution for given problem instance. We follow the frequently taken variational learning approach of using a neural network with parameters 0 to obtain for a given CO problem instance E an associated probability distribution pe( a Ee) over the space of possible solutions 13.. The problem instances E' are assumed to be independently sampled from probability distribution ((E) with Ω supp(g). In this setting the learning problem can be formulated as finding argming Σa q( ((E) Ee Po(a E) E()). Here Eg is shorthand for the sum over all a e 13 In practice both sums are approximated by empirical averages over samples from the corresponding distributions and pa. In CO E(or) is typically characterized by many local minima. Therefore, directly approaching the learning problem described above via gradient descent methods is prone to getting stuck at sub-optimal parameters. As we demonstrate in Fig. (Right) this problem can be alleviated by adding regularization term to the optimization objective Eq. [ which encourages po to R associated to the
avoid a collapse to local minima. The resulting loss function L 111, 1

---

learning problem is given by:
M = M Σ F(0: Fm),
3)
(2)
where M is the number of CO problem instances F. in one batch and | plays the role of regularization coefficient which frequently referred to as the inverse temperature 1/T. The loss for an individual problem instance Em is based the so-called free-energy F(O; B, Em):
F(O; 3, Em) = Σ p (Em + 1/B pe Em))
Em)
(3)
In practice, approximate the expectation value F(0; 3, Fm with empirical sample average F(0; Em) that is based on ng samples ~ (p((oEEm) The term (n) on the right-hand side of Eq. 33 encourages the concentration of po on low-energy solutions. Small values of the regularization parameter E R encourage po(orl Em) to have a large entropy. For a given value of 3 the free-energy in Eq. iis known to be minimized by the Boltzmann distribution associated to E... (see e.g. [MacKay 2003]):
exp-3FFm exp-3Fmm(o)
| 33
(4)
Σ. exp 3 (Emm('))) (o')
Thus by minimizing L the model learns to approximate p (a) Em, 3) for a given Em and B. For 3 - inf the (r(or)Em, has its mass concentrated on the global minima of Em and for 3 - 0 it approaches the uniform distribution [Mézard and Montanari 2009). It is, therefore, to be expected that the minimization of L becomes harder at low temperatures, L.c. as 3 inf, which opens the opportunity to a principled curriculum learning approach (Sec. D. Based on these considerations we reformulate CO problems as the variational problem of approximating P (orl Em 3) in the limit - inf with a variational ansatz Pg that has the variational parameters 6. This problem can be formalized argming limg--uf L(8; ()..
3 Variational Annealing on Graphs
Our method VAG-CO addresses CO as a variational learning problem on graphs. In particular, given set of CO problem instances it learns to approximate the Boltzmann distribution of the corresponding Ising models (Sec. 2 with an autoregressive distribution. To obtain efficient training with this expressive model we apply an annealed entropy regularization. We formulate this learning problem in an RL setting and use PPO to train our model. To alleviate the lengthy sampling process of autoregressive approaches we introduce sub-graph tokenization which allows us to generate multiple solution variables at in a single step without loosing expressive power. We further improve the memory efficiency of our approach by dynamically pruning the graph which represents the CO problem instance.
Autoregressive Solution Generation. In the following specify how represent po(ar FE,) and how to sample it, i.c how to generate solutions σ.
1. Draw problem instance E' (Bi, Jay) ((E) from the data distribution q and construct a graph G = (V, E) based on E. The nodes E V correspond to the spins a; and the set of edges represents the non-zero components of Jıj- The graph G is equipped with node features x; = [B,, f that are associated to its nodes v;, Here 43 is a four dimensional one-hot encoding w which indicates the four possible states (I-IV) of node vi namely whether the corresponding spin o; is set to the value +1 (I) or - (II) or whether it is to be generated in the current step (III) or afterwards (IV). The edges between nodes 1/ and My are associated with the scalar edge features Jij-
2. Order the graph nodes according to the breadth- first search (BFS) algorithm. The i-th node in this ordering denoted o;. Now i and is set to the state (III) and all by are set to (IV).
3. A GNN is used to parameterize Bernoulli distribution (a GNN, (G) from which the value 1 of o, is sampled. The state encoding variables of G are updated by setting
3

---

t associated to the graph is now updated accordingly and 4,41 is set to (III). Now i is
incremented.
4. The previous step is repeated until the values of all are set.
We note that at each step i the graph G depends on the problem instance E and the already generated spins des. Therefore this procedure represents an autoregressive parametrization of a distribution over the space of possible solutions (-1,1}% By denoting the graph G at step G2 we get:
B) I po po E) po(o,)G.).
(S)
This approach is more expressive than MFA and can therefore be expected to be better suited to approximate the typically complex Boltzmann distributions (Eq. 4 encountered in CO. Next specify how we realize a stable training procedure by employing RL methods. The model architecture and hyperparameters are detailed in App. A.10
Reinforcement Learning Setting. To explain how we train our model is convenient to adopt an RL perspective where an episode corresponds to the solution generation procedure described above. For this we consider a Markov Decision Process (MDP) that is given by A, R). Here, the set of possible states and the state at step is denoted by S. At each step s, is represented by the graph G, Given a state s; an action a, that represents the assignments of the spin value a; E (-1, sampled from probability distribution which is parameterised by the policy po Go( After sampling an action the reward R2 (G 99) is computed according to an objective that is based on Eq.
We use the relation F(0; B, F)) Earpp E (,,,,, [3] and define the reward at step i as
1 B, 9)= - |E pog
(6)
where AE, = os [E Joog B,}
With this definition maximising the reward is equivalent to minimizing the free-energy in Eq. We approximate F F(0; 3, E) with the empirical mean of pe(a,la=ss, F.). Finally, the state is changed determin- istically, i.c. the spin is set according to the sampled action a,, The state update corresponds to the updated of tu and f as described in step the solution generation procedure. Our model learns to solve this MDP via the popular PPO algorithm [Schulman et al. 2017 which is described in more
EN R, (3) based on ni samples
detail in App. A.9
Annealing. As discussed in Sec. 2 our training objective is based on Eq. which contains an entropy regularization term. The strength of this regularization is determined by temperature T. At first, the temperature is kept constant at T > 0 for Nwarmop steps. Then, the temperature is slowly adapted for Nanneal steps by following predefined annealing schedule (see App. 1.14) that converges to T 0 as the end of training is reached. This reduction of the temperature throughout training is motivated from a theoretical point of view in Sec.
Subgraph Tokenization. By introducing subgraph tokenization we decrease the number of forward passes per CO problem instance without sacrificing expressivity. Instead of modelling the probability for the two possible values +1 of a single spin, we let the policy represent a probability distribution po(o
over all 2* possible configurations of l consecutive spins in the BFS ordering 2f the solution generation procedure). We represent Pa po(o o E) with a softmax function of a 2 dimensional vector that is output by the GNN (App. modification of step
(App.A..1)) A.10). Subgraph tokenization represents 3 of the solution generation procedure and yields improvement of the performance (Sec.
Dynamic Graph Pruning. We note that once the spin value of o; is sampled its interaction with adjacent spins o, is fixed. Therefore, assuming that we have sampled the first spin 01, its interaction with an adjacent spin o, that is yet to be generated, can be expressed as B(1 - j) where we introduce B11 + f); JI Therefore, we can immediately remove generated spins from the graph and update the node embeddings x, [Bj, (5] of adjacent nodes with B, - B, + B(i + f)j- Reducing the graph size during the solution generation process has the benefit of reducing the memory requirements when processing graphs with a GNN as in step 3 of the solution generation procedure.

---

4 Theoretical Motivation of Annealing
In the following we show that the concept of an annealed entropy regularization in Eq. S can be motivated based the sample complexity of density estimation for Boltzmann distributions.
Consider the problem of approximating the probability density function of Boltzmann distribution (n(s; E, where E R - [o, 1] is a suitably normalized energy function. To solve this task we are given set of samples S = {a1,-..
8m) which are independently drawn from the corresponding target Boltzmann distribution pn((;; E, (8) at a fixed inverse temperature 8 For brevity we will denote the distribution associated to pn( E,33) (3) from now on as p( where E is suppressed since it is a fixed function. The empirical distribution corresponding to S will be denoted as p. Further assume that we can evaluate E(s,) for all s; E E.
In the present context natural feasibility criterion for the estimated distribution is that the expectation value of the energy function for the estimated distribution should be compatible with the corresponding value Bpga (()) of the target distribution. We approach this problem with the maximum entropy principle pJannes 1957]. Informally, this principle prescribes that among all feasible distributions one should prefer the one that yields the highest entropy. Since S contains only finite number of samples we cannot determine pp()) E(s) with arbitrary accuracy. As already stated in a similar form in DDudik et al. 2007 one obtains with Hoeffding's 's inequality that with a probability of 8:
(pp E())))))((() < V 1m(2/0)(22m).
As shown in [Kazama and Tsujii 2003 the resulting maximum entropy optimization problem with the
inequality constraint Eq. 7 is equivalent to the following regularized maximum likelihood problem
over the family of Boltzmann distributions: min, (B) + X3, where (g( = - E, log p(5) is the
cross-entropy between and p(3). Based on closely related result in (Dudik et al. 2007
we obtain
the following bound for the approximation of maximum entropy distributions:
Remark 1. Assume bounded energy function E : RNN- [0. 1] and let 3 minimize Cx(3) 1|31
where x= / /m(22/8)(22m). Then with a probability at least 8:
(xI(p(+))lx(() & |8
2 (2
(8)
See App.AA.1... A.15.2 for further details. According to Eq. the sample complexity of approximating maximum entropy distribution at an inverse temperature 3 is in O(32).
Curriculum learning is machine learning paradigm in which the difficulty of training tasks is gradually increased as training proceeds [Bengio et al., 2009]. In view of Theorem | the entropy annealing in VAG-CO (Sec. 3) can be regarded as principled curriculum learning approach if an increased sample complexity is regarded as an indication of a more difficult learning problem. In supervised learning tasks Wu et al. [ find that the application of curriculum learning results in more resource efficient training but not necessarily in better performance. As we demonstrate in Sec. our entropy annealing does actually yield better performing models.
5 Related Work
As pointed out in (Yehuda et al. 2020 supervised training of CO solvers faces the problem of expensive data generation and that data augmentation cannot circumvent this problem. Consequently, there is a growing interest in RL and unsupervised methods. In the following we focus on methods that attempt to learn how to generate solutions rather than how to improve existing ones.
Unsupervised Learning and Reinforcement Learning. The work of (Bello et al. [017)] proposes to use an actor-critic approach to learn how to solve CO problems. They were the first to show that deep RL is a promising approach for CO. Their method represents an autoregressive distribution over solutions. However, in their setting rewards are sparse since they are only available for complete solutions. In our approach rewards are dense since they are available after every state transition. Another influential work is [Khalil et al. 2017) who were the first to exploit the graph structure of CO problems by learning problem representations with GNNs. Applying GNNs to CO problem has become common choice since then [Cappart et al. 2021 In TTnshooff et al. 2020 the method
5

---

RUN-CSP is introduced. The work of [Karalias and Loukas 2020] (EGN) in which the goal is to minimize loss that can be regarded as a certificate of the existence of CO problem solution whose cost is upper bounded by the loss. By relying on an MFA they can calculate this loss efficiently. Also QQiu et al., 2022] use an MFA to optimize a distribution over the solution space and optimize the corresponding parameters with gradient estimates that are obtained by REINFORCE. The work of [Min et al. 2022] uses the same MFA-basee concept as EGN and focuses on alleviating the oversmoothing problem of GNNs by using modified Graph Convolutional Networks Kipf and Welling 2017]. An approach based on an entry-wise concave continuous relaxation of the discrete CO loss and subsequent rounding procedure for integral solution generation is introduced in [Wang et al. 2022]. Here the concept of [Karalias and Loukas 2020] is generalized to a wider class of problems and rather simple rounding procedure is introduced to generate solutions efficiently. This approach is further extended by [Wang and Li 2023 to a meta-learning setting called Meta-EGN in which network parameters are updated for individual CO problem instances. They also argue that RL based CO methods suffer from unstable training. We call this claim into question by finding no stability issues with RL-based method. Further VAG-CO outperforms EGN and even Meta-EGN despite not updating any parameters on test problem instances. The approach of extending functions on sets continuous supports is taken in Neural Set Function Extensions (NSFE) [Karalias et al.. 2022]. In NSFE the discrete CO loss function is replaced with a convex combination of the discrete loss function values certain integral solutions. These extensions can be regarded as expectation values that can be calculated in closed form without sampling. Whether the lack of sampling noise in NSFE is beneficial for the optimization procedure not obvious.
proposes Erdōs Goes Neural
Variational A Annealing. The concept of using autoregressive models in the variational problem of approximating Boltzmann distributions of Ising models was introduced by (Wu et al., 2019]. They show that variational annealing (VA), i.e. the combination of the variational approach and temperature annealing, is a highly performant method. The work of (Hibat-Allah et al.) [021] compares VA to other ones like Simulated Annealing (SA) (Kirkpatrick et al. 1983] and confirms its strong performance on problems related spin glasses. In [Khandoker et al. 2023) the strong performance of VA compared to SA is confirmed on CO problems. A crucial aspect of these works on VA with respect to ours and the ones in the previous paragraph that they optimize the parameters of their models only for individual problem instances. They do not attempt to learn how to generalize over a family of problems. The work of [Sun et al. 2022] aims at generalizing application of VA in CO by using an MFA. Our experiments indicate that the simplicity of MFA-based methods leads to performance limitations in particular on hard CO problem instances.
6 Experiments
We evaluate VAG-CO on various CO problems that are studied in the recent literature. Additionally, we evaluate VAG-CO on synthetic datasets where solving the corresponding CO problem is known to be particularly hard. Finally, we discuss experiments on the impact of entropy regularization and subgraph tokenization. Our result tables also include inference runtimes. A quantitative runtime comparison is, however, difficult since the runtimes reported from other works and for Gurobi 10.0.0 [Gurobi Optimization, LLC 2023) were obtained with differen setups. See App. 1. for details on the runtime measurements. For Gurobi we report results for various different runtime limits.
Maximum Independent Set. In the following we will compare VAG-CO on the Maximum Inde- pendent Set (MIS) problem to recently published results from [Karalias et al., 2022), where the graph datasets COLLAB, ENZYMES and PROTEINS Morris al. [2200] are used. In the MIS problem the task is to find the largest subset of independent, i.e. unconnected, nodes in a graph. As an optimization objective for VAG-CO we use Eq. 2 with the Ising energy function for MIS (see Tab. D]. Here, the energy function E'q)) consists of two terms E (q) and EB((), that depend on the binary representation of the solution q= a+1 with g e {0, 13 When 4; = 1 the corresponding node is defined to be included in the set and it is excluded otherwise. The first term (Ą(q) is proportional to the number of nodes that are in the set, whereas Eg(q) proportional to the number of independence violations, i.e. the number of connected nodes within the set. By selecting A, B E 14, such that < B we ensure that all minima satisfy En = since this case excluding a violating node from the set always reduces the energy. In our experiments, we choose A -= 110 and B 1.1.
We follow the procedure of Karalias et al. 2022] and use a 2...6/0.///0 train/val/test split on the aforementioned datasets and use only the first 1000 graphs on the COLLAB dataset. The results

---

Problem Type Ising formulation: min, E(g) where E(q) = Ex(g)= + En(g)
MVC
E(q) - A E.
q: + B. E Denneel (1 q)) (1 (1-99) qs)
MIS
k(g) = - A Mix 9 + B Disie 9:
MaxCI
E(q)- = - - L. qi + B Lasne 9: 4j
MaxCut
E(o) = -- -Eisee 1-g. where 2q
Table 1: Ising formulations of the MVC, MIS, MaxCI and MaxCut problem (Lucas 2014)) The term that includes the constant A corresponds to (g) and the term that includes the constant B corresponds to Eg()).
ENZYMES
PROTEINS
IMDB BIDB-BINARY
COLLAB
MUTAG
Evalustion metric A.1
AR (s/graph)
AR" (s/graph)
ARR (s/graph) )
AR" (s/grphh) 1
AR" (s/graph
(e)
0.821 0.124 (N/A)
0.903 0.114 (N/A)
1: 0.310 (N/A)
0.886 0.198 (N/A)
0,589
(N/A)
NSFE
0.775 0.115 (N/A)
0.729 0.205 (N/A)
0.679 1: 287 (N/A)
0.592 0.253 (N/A)
0.854 1: O. 152 (N:A)
REINFORCE (r)
0.751 1 0.301 (N/A)
0.725 :: 0.285 (N/A)
0.881 1: 0.240 (N/A)
1.000 (N/A)
0.781 1: 0.316 (N/A)
Straight through
0.725 0.268 (N/A)
4: (N/A)
0,253 (N:A)
0,856 0.221 (N/A)
0,065 0.302 (N:A)
DB-Groedy
0.5810
0.0000 (0.008)
0.9848 + 0.0005 (0.008)
1.000 (0.007)
0,5068 0.0001 (0.033)
0.384 + 0.061 (0.006)
MFA:
0.9675 0.001 (0.036)
0.9 0.9683 :: 0. 0007 (0.037)
1.000 (0.022)
0. 0.9065 0.0004 (0.054)
1.000 (0.019)
MFA-Amnall CE
0.9880
0.0014 (0.036)
0,9681 0. 0007 (0.037)
1,000 (0.022)
0.9991 0.0004 (0.054)
1.000 (0.019)
VAG-CO (vurs)
0.9960
+ 0.0007 (0.026)
0.9977 + 0.0005 (0.032)
1.000 (0.017)
0.9088 0.0002 (0.064)
1.000 (0.016)
Gurobi Moss 0.001)
1000 (0:0011)
0.8866 + 0.0001 (0.001)
1,000 (0.0003)
0.997 0.001 (0.002)
1,000 (0.0001)
Gurobi (h..
(0.001))
1.000 (0.001)
1.000 (0.0003)
1.000 (0.002)
1.000 (0.0001)
Evalustion motric A.1
Alt (s/srapla)
Aît (s/arapha >
(s/graph
Aft (s/graph)
AR is/raph
MFA
0,9635 0.0014 (0.0029)
0,9705 0.0012 (0 0025)
0.990 +: 0.002 (0. 00288
0.960 0.007 (0.0024)
0.970 0.003 (0. 0028)
MFA-Ammeal
0 + 0.001 (0.0029)
0.073 = 0.002 (0.0025)
0.984 + 0.001 (0.0028)
0.917 = 0.045 (0.0024)
0.976 + 0.002 (0.0028)
VAG-CO (purs)
0,9863 0.0007 (0.026) 0.9908 o. 0.0005 (0.032)
0.905
(0.017)
0.001 (0.054) 0.994 0.002 (0.015)
Table 2: Results on the MIS Problem (see Sec. D. We show test set results on the best approximation ratio A and the average approximation ratio AR across different methods and datasets. that are closer to The time for each algorithm is reported in round brackets as seconds per graph.
Values
better. (r) indicates that these results reported in [araliaa et al. 2022
are shown in Tab. 2 where the test set average of the best approximation ratio AR out of ns 8 sampled solutions per graph is reported (see App. A.I). This metric was originally proposed in [Karalias and Loukas 2020). We also report results of our own implementation of Maseb
method that is trained with REINFORCE. This method is used with (nFFAA) and without
(MFA) annealing (App. A.5 p. As in Wang and Li [2023] and Karalias and Loukas (2020) we use the conditional expectation procedure (CE, App./App....... A.5.5) to sample solutions and report the corresponding results as (MFA: CE) and (MFA-Anneal: CE). We also add results obtained by the Degree Based
Greedy (DB- (DB-Greedy) algorithm [Wormald 1993) that is proposed by (Angelini and Terci-Tirs 2023) as a baseline for machine learning algorithms on MIS. For VAG-CO we greedily sample different states for given problem instance by using different BSF orderings of the corresponding graph nodes. Our results show that VAG-CO significantly outperforms all competitors including the MFA on ENZYMES and PROTEINS. On IMDB-BINARY and MUTAG the MFA-based approaches and VAG-CO outperform all other machine learning methods and achieve an optimal AR* On COLLAB MFA-Anneal and VAG-CO exhibit the best results with insignificantly better results for MFA- Anneal. We also report results based on the test set average of the approximation ratio with ns = 30 AR (App. A. This metric shows the performance of the learned probability distribution for each model, when no post processing applied. Our results show that VAG-CO always achieves a significantly better performance in terms of AR than based approaches.
2023 where the Minimum Vertex Cover (MVC) problem is solved on the TWITTER [Leskovec and Krevl 2014]. COLLAB and BOOB-BINARYII [Morris et al. [2200] graph datasets. The MVC problem is the task of finding the smallest subset of vertices such that each edge has at least one node in the subset. As for the MIS problem we formulate this CO problem in terms using the Ising energy function that defined in Tab. [ and set A 1.0 and B 1.1. We follow the procedure of Wang and Li 2023 and use 0.7 0.7/O.I, / train/val/test split on these datasets. Their method Meta-EGr uses meta
Minimum Vertex Cover. We compare VAG-CO to results from Wang and Li)
7

---

learning in combination with the EGN method. They apply CE to generate solutions and report the best out of Ss = 8 sampled solutions. Additionally, they report fine tuning results denoted by (f-t), where they apply one step of fine tuning individually for each CO problem instance in the test dataset. Our results in Tab. show that VAG-CO improves upon the results reported in Wang and Li| (2023) across all datasets. We also show that the MFA with and without annealing performs remarkably well on the TWITTER and COLLAB dataset and even outperforms VAG-CO on the COLLAB dataset. The results in terms of AR show that the strong performance of the MFA can be attributed to the post-processing method CE rather than to the learned MFA- -rased model.
Maximum Clique. The Maximum Clique (MaxCI) Problem is the problem of finding the largest set of
TWITTER
COLLAD
IMDB-BINARY
Evaluztion metric A.1
AR" (s/gaph))
ARC (s/graph)
ARe (s/graph)
EGGN
LAMH 0.02 (0.29)
1,013 0,022 (0.15)
LOUD [0.06]
DGN F-1
1,028 a 0.021 (0.80)
1,008 1. 0.015 (0.38)
1.000 (0.32)
nodes, in which each node is connected to all other nodes that are contained in the set. The MaxCl prob-
RUN-CSP
1.100 1 0.435 (0.
1,208 1: 0.198 (0.10)
1188 1: 0.178 (0.08)
Meta-EGN o1
1.119 :: 0.000.00.0200
1,003 1. 0.010 (0.15)
1.000 (0.08)
Mete-EGN (+)
L017 0.017 (0.80)
1,002 1: 0.010 (0.35)
1.000 (0.32))
Groady (1)
0.014200014:.000
1,200 + 0. 0.105 (1.70)
1,180 + 0,077 (0.02)
MFA: CB
I 10004 0.0000 (0.08)
100021 1 0.00003 (0.06)
1.000 (0.02)
lem is equivalent to solv-
Amnel CE
1 1,0011
:: 0.0001 (0.08)
1,00255 00004 (0:05)
1.000 (0.02)
VAG-C0 lours)
1,004 0.00006 (0.15)
1,0017 1: 0,0002 (0.
1.000 (0.016)
ing the MIS problem on the complementary graph. The complementary graph can be obtained by con- necting nodes that are not connected in the original graph and by disconnecting nodes that are connected. We evaluate our method on the ENZYMES and IMDB- BINARY datasets with the same train/val/test splits that are used in the MIS problem (see Sec. 6). Results are shown in Tab. where we compare to results that are reported in [Karalias et al.
Gurobi (fa.. 0.01)
1.884 :: 0:001 (0.008)
1,002 :: 0,0001 (0.0022
1.000 (000005)
Gurobi
1:0055 0.00001 (0.125)
1000 (0.000))
1.000 (QUOOIS))
Gurobi Urane 0.2
1.0001
0.0001 (0.028)
1.000 (0.004)
1.000 (0.0000)
Evalustion metricA. A.1
AR (s/mmph))
AR (s/graph)
AiR (s (s/grapha
MFA:
10106 0.0002 (0.004)
10010 0.00001 (0.000)
1.014
0.0004 (0,005)
MIA-Anned Anneal:
1.0000 0.0001 (0.004)
1.0048
O. 0.0002 (0.003)
0,0002 (IL/K33) 1.001
1.0011
VAG BAG-CO (ours)
1.0079 0.0003 (0.15)
1,0056 1. 0.0002 (0.16)
0,0006 (0.016)
Table 3: Results the MVC problem (see Sec. D.
We show test set results on the best approximation ratio AR* and the average approxi- mation ratio AR across different methods and datasets. Values that are closer to one are better. (r) indicates that these results are reported in (Wang and Li 2023).
2022 and to our own implementation of MFA and MFA-Anneal. Here, we see that VAG-CO is always among the best performing method in therms of AR* and always the single best method in therms of AR.
Maximum Cut. The Maximum Cut (MaxCut) problem is the problem of finding two sets of nodes so that the number of edges between these two sets is as high as possible. We evaluate our method on BA graphs [Barabási and Albert 1999]. where we follow (Zhang et al., 2023] and generate 4000/500/500 train/val/test graphs between the size of 200-300 nodes and set the generation parameter m to 4.
Results on this dataset are shown in Tab. where we compare our method MaxCut Values
(MCut E(EJYe 2
I 1-a,a; that are reported in (Zhang et al. 2023). Our method achieves the best
MCut results among all learned methods.
1.000
1.035
10F
0.975
1,030
MiA CE L--33
Gandl 122
MJ 7-- 10
0.050
0.925
1.020
508-45
wws.wwarn 133
Gambl
VAG-CD with ST 7g-
VAG CD m/
VAG-CO 5
M/ -6--00-
0.900
1.015
2 303
0.875
1.030
0.850
0.825
Oodi
1.005
0.800
5.0
10.0
12.5 15.0 17.5 20.0
0.3
0.4
number of gradient
Figure 1: Left: Comparison of AR for the MIS problem on the RRG-100 dataset with different degrees d. Middle: Results for the MVC problem in terms of AR* on the RB-200 dataset at different generation parameters p. Left, Middle: Values closer to one are better. The runtime t is given in s/graph and fmax the time limit for Gurobi. Right: Ablation on annealing and subgraph tokenization (ST). The Evol on the RB-100 MVC is plotted over the number of gradient update steps. Lower values are better.
8

---

ENZYMES
IMDB BINARY
BA 200 300
Evaluation metno A.1
AR (s/ graph)
ARR (s/goaph))
MCai (s/graph)
EGN (1)
0.883 = 0.156 (N:A)
0338 0.176 (N/A)
NSFE (
0.933 48 (N/A)
0.961 0. 0.143 (N/A)
Gurobi (r)
732.47 (1.57)
REINFORCE (1)
0.751 # 0.301 (N:A)
0.881 0.240 (N/A)
SDP (r)
700.36 (4.29)
Straight-Bnoighhh (1)
0.725 +: 0.268 (N/A)
0.917 253 (N/A)
Greedy (O)
688.31 (0.026)
DB Greedy
0.9685 1. 0.0035 (0.025)
0,9873 0.0008 (0.01)
BON (f)
693.45 (0.0892
MEA; CB
0,9765 0.001 (0.04)
0.999 0.009 (0.125)
MFA-Anaeal: CE
0.9921 = 0.0017 (0.04)
0.999
0.0001 (0.025)
Anneal (f)
696.73 (0.09)
VAG-CO (ours)
0,987 0.004
0.9981 0,0013 (0,917)
GFlowNet (f)
704,30 (0.35)
Garobi Crase 0.01)
0.983 + 0.002 (0.004)
0.596 + 0.0002 (00001)
VAG-CO (ours)
722.31 (0.17)
Gurbi Unaee
(0.005)
1,000
Gurobi (fmaxe 0.1)
731.34 (0.1)
Evaluation metric
AR in (e/graphi
AR
(e/groph))
Gurobi (Omex 1.)
731.79 (1.)
MEVA-Anneal
MFA
0.764 0.013 (0.001)
0,806 0.001 (b:001)
0.983 + 0.009 (0.000)
0 + (0.001)
Gurobi (fmax 2.)
732.00 (2.)
VAG-CD
0.943 0.012 (0.029)
0.988 0.0013 (0.017)
Gurobi (fmax = 10.)
732.01 (10.)
Table 4: Left: Results on the Maximum Clique Problem (see Tab. 1. We show test set results on the best approximation ratio AR and the average approximation ratio AR across different methods and datasets. Values that are closer to one are better. (r) indicates that these are results as reported in [arallaa et al., 2022]. The time for each algorithm is reported in round brackets as seconds per graph. Right: Results on the Maximum Cut Problem (see Tab. [ on BA 200 - 300 graphs. We show
test set results on the average maximum cut value MCut across different methods. Values that are larger are better. (r) indicates that these are results as reported in [hhang et al. 2023
Evaluation on Synthetic Problems for MIS. On many of the benchmarks above nearly optimal results are obtained. Therefore, we conduct additional experiments on synthetic graph datasets that yield hard CO problems.
For graphs with a degree d larger than 16 the MIS problem on random regular graphs (RRGs) is known to be particularly hard [Barbier et al. 2013 Therefore, we generate RRGs with an average size of 100 nodes (RRG-100), where we sample 4100 RRGs between the size of 80 to 120 with different degrees d E [3, 7, 10, 20]. For training, validation and testing we use 3000/ 1000//0///// 00/1000 graphs. Results for the MIS problem on the RRG- 100 dataset are shown in Fig. Fig. (Left), where we plot the test set average of the best relative eπror cet out of ng 8 sampled solutions per problem instance. Here the relative error is calculated with respect to exact solutions that are obtained with the commercial solver Gurobi. Error bars indicate the standard error over graphs with the corresponding node degree. As expected we find that for all methods rrel increases for larger degrees. MFA and MFA-Anneal outperform the DB- OB-Greddy method. Furthermore, VAG-CO is the best performing method for all graph degrees. We also show Gurobi performance for each d at different time limits mmax- On this dataset VAG-CO outperforms Gurobi with the given compute budged (see App. A.12) by a high margin on hard instances.
Evaluation Synthetic Problems for MVC. We also conduct experiments on graphs that are generated by the so-called RB method [Xu al., 2005]. This method allows the generation of graphs that are known to yield hard MVC problem instances [Wang and Li 2023). The RB model has three distinct generation parameters: n, k', and p. Adjust- ing the values of n and k allows us to control the expected size of the generated graphs, while serves as a parameter to regulate the hardness of the graph instances. Specifi- cally, when is close to one, the generated graphs tend to be easier, whereas reducing p leads to increased difficulty. To ensure diversity in the RB dataset, we generate graphs with varying levels of hardness by randomly sampling p from the range of 0.25 to 1.0. For our experiments, we utilize the RB-200 dataset, consisting of RB graphs with an average size of 200 nodes with a train/val/tess split of 2000/500/ 5000/0000 graphs. Following [Wang and Li, 2023) each graph in this dataset is generated by randomly select- ing values for n from the range of 20 to 25, and k'' from the range of 9 to 10. Figure| (Middle) shows the corresponding results in terms of AR. in dependence
1.06
0.4
1.04
S0reeph
1.03
0.2
1.02
2
10
Subgraph Lokenication k
Figure 2: Study on the advantages of larger ST configuration sizes k for VAG- CO on the RB-100 MVC dataset. On the left y-axis we plot AR* (lower val- ues are better) for VAG-CO models that are trained for values of k. The left y- axis shows the inference time per graph.

---

of p. We find that the DB- DB-Greedy algorithm outperforms both MFA approaches. The MFA with annealing tends to perform worse than the MFA without annealing on this dataset. Here, we also add our own implementation of EGN and EGN-Anneal (Sun et al. 2022 and find that EGN performs similarly to MFA and EGN-Anneal achieves a similar performance as DB-Grest OB-Greedy. VAG-CO out- performs all other methods for all values of on this dataset. We also show Gurobi performance for each p at different time limits fmax. Here, we see that at high p values Gurobi achieves close to optimal results within very small time budged and is the best performing method. Whereas for low p values Gurobi performance drops by a high margin and VAG-CO achieves similar results within a comparable time.
Ablation on Annealing. In the following we investigate whether the annealed entropy regularization in Eq. 3 is indeed beneficial. These experiment are conducted for the MVC problem on RB-100
graphs that are generated with e [9. 15]. k' [8, 11] and E [0.25. 1]. We compare VAG-CO
2, 2,5 10 3 03. The run with To = 0 is equivalent to VAG-CO without annealing (App. A.I1.. Figure (right) shows the average relative error érel for ris = 30 sampled solutions per problem instance on the test set over the number of gradient update steps. We find that the run without annealing starts to converge rather quickly but at a comparably bad relative energy. In contrast to that, the runs with annealing (To (5 10 2,5 10 10=3)) keep improving over more steps and achieve better final performance.
test set learning curves with three different initial temperatures To E {5 100211020000
Ablation on Subgraph Tokenization. Next we study the impact of subgraph tokenization (see Sec. 3 on the performance of VAG-CO. The same problem setting as in the annealing ablation above used. Figure | (right) compares VAG-CO with the default k = subgraph tokenization (w/ ST) to VAG-CO without subgraph tokenization (w/o ST). For a fair comparison the numbers of gradient update steps and the annealing schedules are set to be equal. We find that with subgraph tokenization we achieve a considerably better érel. These results underpin that subgraph tokenization does indeed yield an improved efficiency in terms of gradient update steps.
We additionally provide a more detailed investigation on ST. where we compare VAG-CO with four different values of & E {1, 2, 5, 10} on the RB-100 MVC dataset. In our experiments we keep the number of gradient update steps constant and iteratively tune for the best learning rate and the best initial temperature in each setting of k. Results are shown in Fig. [ where AR and the average time per graph are plotted over k. As k increases the performance in terms of AR* does improve and the inference time reduced considerably. This experiment underscores the performance and scalability improvement due to ST.
Limitations
While MFA-based approaches typically generate solution CO problem in single forward pass VAG-CO requires a number of forward passes that is proportional the problem size N. However, methods with solution generation in single forward pass frequently rely on post-processing proce- dures like CE that have a time complexity which is also linear in N. In our approach the nodes in the graph are always processed in order that is determined by BFS. Therefore studying the impact of other node orderings or making the algorithm invariant to such orderings might be an interesting future research direction. For the annealing schedule report comparisons for VAG-CO with three different initial temperatures but comprehensive study on the optimization of the annealing schedule is left to future investigations.
8 Conclusion
Our results show that learning solve hard CO problems requires sufficiently expressive approaches like autoregressive models. With VAG-CO we introduce method that enables stable and efficient training of such a model through an annealed entropy regularization. We provide a theoretical motivation for this regularization method and demonstrate its practical benefit in experiments. In addition, we demonstrate that by systematically grouping solution variables VAG-CO achieves an improved performance with respect to a naive autoregressive model. Importantly VAG-CO outperforms recent approaches on numerous benchmarks and exhibits superior performance on synthetic CO problems that are designed to be hard.
10

---

Acknowledgements
We thank Günter Klambauer for useful discussions and comments, The ELLIS Unit Linz, the LIT AI Lab, the Institute for Machine Learning, are supported by the Federal State Upper Austria. We thank the projects AI-MOTION (LIT-2018-66 YOU-212), DeepFlood (LIT-2019-33 2--00-99--70 YOU-213), Medical Cognitive Computing Center (MC3), INCONTROL-RL (FFG- 88 1064). PRIMAL (FFG-873979) S3AI (FFG-872172), DL for GranularFlow (FFG-871302) EPILEPSIA (FFG-892171), AIRI FG 9-N (FWF-36284, FWF-36235), AI4GrrennHeaiingGinssFFFCC 899943), INTEGRATE: (FFG-892418), ELISE (42000-ICT-2019-3 ID: 951847), Stars4Waters
WWWT:ONN-C.--2211-LLLMAA-AA We thank Au Deep Learning Center, TGW LOGISTICS GROUP GMBH, Silicon Austria Labs (SAL), FILL Gesellschaft mbH, Anyline GmbH, Google, ZF Friedrichshafen AG, Robert Bosch GmbH, UCB Biopharma SRL, Merck Healthcare KGaA, Verbund AG, GLS (Univ. Waterloo) Software Competence Center Hagenberg GmbH, TÜV Austria, Frauscher Sensonic, TRUMPF and the NVIDIA Corporation.

---

References
Ralph Abboud, ismail ilkan Ceylan, Martin Grohe, and Thomas Lukasiewicz. The surprising power of graph ncural networks with random node initialization. In Zhi-Hua Zhou, editor, Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI 2021, Virtual Event / Montreal, Canada. 19-27 August 2021, p 118.2218. ijcai.org, 2021. doi: 10. 201222991 URL https://doi org/ 10. 21963/29943/41 h/t/////jaai../ 2021/291
Abien Fred Agarap. Deep learning using rectified linear units (relu). CoRR, abs/1803.088375 2018. URL. httt://arxivg org/ http///aaa//./e...e// 08375
Sungsoo Ahn, Younggyo Sco, and Jinwoo Shin. Learning what to defer for maximum independent sets. In Proceedings the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event, volume 119 of Proceedings of Machine Learning Research, pages 134-144. PMLR, 2020. URL
http://proceeddnnge mlr. www...pressvvviv arees//v11/ahhh0a.a
Maria Chiara Angelini and Federico Ricci- Tersengii Modern graph neural networks do worse than classical greedy algorithms solving combinatorial optimization problems like maximum independent set. Nat. Mac. Intell., 5(1):29-31, 2023. doi: 10. 00010088/62255--02-0005999 URL https://doot otg/2o..g111.. 1038/
-42-56-022-000
Lei Jimmy Ba, Jamie Ryan Kiros, and Geoffre F. Hinton. Layer normalization. CoRR, abs/1607.064550, 2016. URL httpt://///aiii org/ 1607 1607006450 06450
Albert-László Barabási and Réka Albert. Emergence of scaling in random networks. Science, 286(5439): 509-512, 1999. doi: 10.1 122/sccincce286.5555559 1126/ 286 52395499555 509
URL hittps://ww/ httpps:/www..acience. science. www.eeeneeeereeeeeeeeeeee wrtp/oi/a//////10 abs/10
Jean Barbier, Florent Krzakala, Lenka Zdeborová, and Pan Zhang. The hard-core model on random graphs revisited. Journal of Physics: Conference Series, 473(1):012021, dec 2013. doi: 10. 10.1088/1742-659661 /317/////66///////////
012021. (RR.sstts:://dt http:://xxxddddd..ooo/ www.oo.......ooooo 111.000//1144--006////////1/00/0
Peter W. Battaglia, Jessica B. Hamrick, Victor Bapst, Alvaro Sanzhez-Gonzaler,, Vinícius Flores Zambaldi, Mateusz Malinowski. Andrea Tacchetti, David Raposo, Adam Santoro. Ryan Faulkner, Çaglar Gülçehre. H. Francis Song, Andrew J. Ballard, Justin Gilmer, George E. Dahl, Ashish Vaswani, Kelsey R. Allen, Charles Nash, Victoria Langston, Chris Dyer, Nicolas Heess. Daan Wierstra, Pushmeet Kohli, Mathhw M. Botvinick, Oriol Vinyals, Yujia Li, and Razvan Pascanu. Relational inductive biases, deep lcarning, and graph networks. CORR, abs/ 03//1006.1161, 2018 http:///aaax http:///aaxivv..gaabbee/// http//880600.1..00
Irwan Bello, Hieu Pham, Quoe V. Le, Mohammad Norouzi, and Samy Bengio. Neural combinatorial optimization with reinforcement learning. In 5th International Conference on Learning Representations, ICLR 2017. Toulon, France, April 24-26, 2017 Workshop Track Proceedings. OpenReview.net, 2017. URL. https
//openreviev.
Yoshua Bengio, Jérôme Louradour, Ronan Collobert, and Jason Weston. Curriculum learning. In Proceedings of the 26th Annual International Conference on Machine Learning. ICML '09, page 41 48, New York. NY, USA, 2009. Association for Computing Machinery. ISBN 678160558111.. doi: 10. 15531555533744553888
URI htt.t::/p///oii orgg 10. 1145/ 005555555337.... 1553380
Yoshua Bengio, Andrea Lodi, and Antoine Prouvost. Machine leaming for combinatorial optimization: A methodological tour d' hororizon. Eur. J. Oper. Res., 290(2):405-42 2021. doi: 200100606.
URL /atts://s:///oi. http://ooo.org//// org/10. wttg//....///..............//..
2020. 00.00...... 066
Maximilian Böther, Otto Kißig, Martin Taraz, Sarel Cohen, Karen Seidel, and Tobias Friedrich. What's wrong with deep learning in tree search for combinatorial optimization. In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net, 2022. URL ? 7id Y=mkkHzzgq7111
httpps:///penreiee..eeee
Quentin Cappart, Didier Chételat, Elias B. Khalil, Andrea Lodi, Christopher Morris, and Petar Velickovic. Combinatorial optimization and reasoning with graph neural networks. In Zhi-Hua Zhou, editor, Proceedings of the Thirtieth International Joint Conference Artificial Intelligence, IJCAI 2021, Viriual Event / Montreal. Canada, 19-27 Augusi 2021, pages 4358-4355. ijeai.org, 2021. doi: 0..24963/ieaii222215553 URL URt/ttpps //doi //doi.org/10 org/10. 10. 24963/ ijcaiiicaaaaaa 2021/595
Corinna Cortes, Vitaly Kuznetsov, Mehryar Mohri, and Umar Syed. Structural maxent models. In Francis Bach and David M. Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, ICML 2015, Lille, France, 6-11 July 2015, volume 37 of JMLR Workshop and Conference Proceedings, pages 391-399, 399. JMLR.org, 2015. URL Uttprptpp// mlr wwwpppreesr// v37////////// ctr//oortesss/ htmll
12

---

Miroslav Dudik, Steven J. Phillips, and Robert E. Schapire. Maximum entropy density estimation with generalized regularization and application species distribution modeling. Mach. Learn. Res. 8: 1217-1260, 2007. doi: (1.555//3114998.1144554 URL https: http://ddlaaa..././../ 5655///11..555111449999999 1314540
Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2023, URL https:://ww..gtuue www.wggrbb..... com
Mohamed libat-Allah Estelle M. Inack, Roeland Wiersema, Roger G. Melko, and Juan Carrasquilla. Variational neural annealing. Nat. Mach. Intell., 3(11):952-99) 2611, 2021. doi: 10. 101000/422560002-00000000
URL
htttps://oi. org 10. 1038/ /42.//42256--o 021 - 00401-3
Tommi S. Jaakkola and Michael I. Jordan. Improving the mean field approximation via the use of mixture distrihutions. In Michael Jordan, editor, Learning Graphical Madels, volume 89 of NATO ASI Series, pages 163-173. Springer Netherlands, 1998. doi: 111007/988-94--1-5550--- org/ 10. 10070...107978818 0000007977894--11100 04-01115014-9
URL https://doi
E. T. Jaynes. Information theory and statistical mechanics. Phys. Rev, 106:620-633 630. May 1957. doi: 10.1 PhysRev.106.620. URL aattpa://link. aps org/ http://oo//...1///ttttee..eee 1103/PhysRev....mm 106. 620
Nikolaos Karalias and Andreas Loukas. Erdos goes neural: an unsupervised learning framework for combinatorial optimization graphs. In Hugo Larochelle, Mare' Aurelio Ranzato, Raia Hadsell, Maria-Florina Balcan, and Hisuan-Tien Lin, editors, Advances Neural Information Processing Sys- tems 33: Annual Conference on Neural Information Processing Systems 2020. NeurIPS 2020, De- cember 6-12, 2020, virtual, 2020. URL htttps://proceeddnng wwurips.cc/papeeeeeee wttp:cccpppppppeeeeeeeeeeee /2t00//as// Abstract.htm
Nikolaos Karalias, Joshua Robinson, Andreas Loukas, and Stefanic Jegelka. Neural set function extensions: Learning with discrete functions high dimensions. NewrIPS, 2022. URL http://pppppppr inaappers.iipppp. mips.ipp..... cc/
html
ichi Kazama and 'cch Tsujii. Evaluation and extension of maximum entropy models with incquality constraints. In Proceedings of the Conference Empirical Methods in Naturul Language Processing, EMNLP 2003, Sapporo, Japan, July 11-12, 2003, 2003. UtL.lttts:://clleaaa hctp///laalnttoogg.....oo org/W03-1018/
Elias B. Khalil, Hanjun Dai, Yuyu Zhang, Bistra Dilkina, and Le Song. Learning combinatorial optimization algorithms over graphs. Isabelle Guyon, Ulrike von Luxhurg, Samy Bengio, Hanna M. Wallach, Rob Fergus, S. V. N. Vishwanathan, and Roman Garnett, editors, Advances in Neural Information Processing Systems 30: Annual Conference Neurul Information Processing Systems 2017, December 4-9, 2017 Abstract.
Long Beach, CA. USA, pages 6348-6344 6358. 2017.
neurips ht/p:c/ppppaaee. hac/pappr/20111
Shoummo Ahsan Khandoker, Jawaril Munshad Abedin, and Mohamed Hibat- Allah. Supplementing recurrent neural networks with anncaling to solve combinatorial optimization problems, Mach. Learn. Sci. Technol. 4(1); 15026, 2023. doi: 10. 10118866222222 25533b8895.. URL. http:::::/ii..rr////////// 10tp://1...188/2221110 21t///222222132/2222322322 acb895
Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In International Conference Learning Representations, 2017. httptppp///ggeeeeieee.eeeeet///t Yd=SJU4
S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by simulated annealing. Science, 220(4598): 671-680, 1983. doi: 10.1 2tt.science.22.44598... URL http::tttt/....ciiiittttae wwwpwwwwwciinncee...... www.annc..ooooo.o...oeeeoooooooe org/doi/abs/10 1126/science 200.22224550 200.4598.677
Wouter Kool, Herke van Hoof, and Max Welling. Attention, leam to solve routing problems! In 7th International Conference on Leurning Representations, ICLR 2019, New Orlewns, LA, USA, May 6-9, 2019. OpenReview.net, 2019. h
Jure Leskovec and Andrej Krevi, SNAP Datasets: Stanford large network dataset collection. http://snap stanf franford.r www.oorddddddddaaaaeeeeeeeeeaeae edu/data] June 2014.
Fenglin Liu, Xuancheng Ren, Zhiyuan Zhang, Sun, and Yuexian Zou. Rethinking skip connection with layer normalization transformers and resnets. CoRR, abs/2205500205, 2021. URL http:://////a
2105.07205
Ilya Loshchilov and Frank Hutter. SGDR: Stochastic gradient descent with warm restarts. In International Conference Learning Representations, 2017. URL httpp:///penrrewier http://s///.pppreeeee/eee//////a
13

---

Andrew Lucas. Ising formulations of many up problems. Frontiers in Physics. 2, 2014. ISSN 2296-424X. doi: 00.338///hhy20144.0000 URL htttttt/://we frontiersin www.:fttniiiiiooorrrrrreeeeeee.e wwwporraigaaae/iirrooooooor wwtp:////.................. 2014 00005
David J. C. MacKay. Information theory, inference, and learning algorithms. Cambridge University Press, 2003.
ISBN 978-0-521-64298-9
Nina Mazyavkina, Sergey Sviridov, Sergei Ivanov, and Evgeny Burnaev. Reinforcement learning for combinato- rial optimization: A survey, Comput. Oper. Res.. 134:105400, 2021. doi: 10. 100161i.co:2021110554000 URL
https://doi. org 10. 1016/j
2021. 105400
Song Mei, Theodor Misiakiewicz, and Andrea Montanari. Mean-field theory of two-layers neural networks: dimension-free bounds and kernel limit. Alina Beygelzimer and Daniel Hsu, editors, Conference on Learn- ing Theory, COLT 2019, 25-28 June 2019, Phoenix, AZ. USA, volume 99 of Proceedings of Machine Learning Research, pages 2388 2464. PMLR. 2019. http:tttt:///ceeeeeeeeeeee
mlr wwwppeewseeeepeepe http:/9999/mmmm9m// ht///me19... html
Yimeng Min, Frederik Weakel, Michael Perlmutter, and Guy Wolf. Can hybrid geometric scattering networks help solve the maximum clique problem? In NewrlPS, 2022. URL httt:://papers mips.ipppp... www.nicpppppppeeeeee.eep caper.
http:/0000000000000001001000aaa0 htw.:/4407774----staaa--------e- Conferrre html
Christopher Morris, Nils Kriege, Franka Bause, Kristian Kersting, Petra Mutzel, and Marion Neumann.
Tudataset: A collection of benchmark datasets for learning with graphs. CoRR, abs/2007.08663, 2020. URL
https://arxiv org/ abs/ 2007 08663
Mare Mézard and Andrea Montanari. Information, Physics, and Computation. Oxford University Press, 01 doi: 10. 10000///amm00o007777700777777000 10. 10t11093/accac /0o0///cpoof..... : 001.0001
2009. ISBN 9780198570837,
URL httptttttppp.o...ooo org/
Christos Papadimitriou and Kenneth Steiglitz. Combinatorial optimization: algorithms and complexity.
Courier Corporation, 1998.
G Parisi. Statistical field theory, vol. 66 in. Frontiers in Physics, page 343, 1988.
Ruizhong Qiu, Zhiqing Sun, and Yiming Yang. DIMES; A differentiable meta solver for combinatorial optimization problems. In NewrlPS, 2022. http:///pp/paape nips cc/ sasper_fflle/eeeeeeeeeeeeeeseeee 2022/
html
Prabhakar Raghavan. Probabilistic construction of deterministic algorithms: approximating packing integer programs. Journal of Computer and System Sciences, 3700)1300114, 1988.
John Schulman, Philipp Moritz, Sergey Levine, Michael I. Jordan, and Pieter Abbeel. High- High-dimensional continuous control using generalized advantage estimation. In Yoshua Bengio and Yann LeCun, editors, 4th International Conference Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings, 2016. URL http://tt:///rt: org /ba//506...... .22488
John Schulman, Filip Wolski, Prafulla Dhariwal, Alce Radford, and Oleg Klimov. Proximal policy optimization algorithms, CoRR, abs/1107.0632 00337777777..... 2017. URL http://ttaxxii>> org/ awspaaa/.///////// 06347
Haoran Sun, Etash Kumar Guha, and Hanjun Dai. Annealed training for combinatorial optimization on graphs. In OPT 2022: Optimization for Machine Learning (NeurIPS 2022 Workshop), 2022. URL https //openreview. oet. TkU
Jan Tönshoff. Martin Ritzert, Hinrikus Wolf, and Martin Grobe. Graph neural networks for maximum constraint satisfaction. Frontiers Artif. Intell., 3:580607, 2020. doi: 15338//fri.2220.5550 URL httt:::/t///o
org/10. 2...0389/fffe frai. 2020. 580607
Martin J. Wainwright and Michael I. Jordan. Graphical models, exponential families, and variational inference. Found. Trends Mach. Learn.. 1(1-2):1-3000 305, 2008. doi: 10. 10.1561/220000001, https:////ddo org/ 10 1561/2200000001
Haoyu Wang. Nan Wu, Hang Yang. Cong Hao, and Pan Li. seare for combinatorial optimization with principled objective relaxation. In NewrIPS, 2022. URL http:://pppprs..ttteea mips httppp.cccpppp...pee fomer_fflee files/ paper/2022, hate////hah//ae cbc wbc.aad206fffaafae aww.aaabbbbbbbee00errrrer o0 orw.raaaseseeeeeeaeaaaaaaaaaaa-a Abstract- Conferenccccceeeeeeoeoe html
Haoyu Peter Wang and Pan Li. Unsupervised learning for combinatorial optimization needs meta learning. In
The Eleventh International Conference on Learning Representations, 2023, URL https:///oenreviee
met/fooum?id- EEW.CdE22EEECEBEB
14

---

Nicholas C Wormald. Differential equations for random processes and random graphs. The annals of applied probability, pages 1217-1235, 1995.
Dian Wu, Lei Wang, and Pan Zhang. Solving statistical mechanics using variational autoregressive networks. Phys. Rew. Lett.. 122:080602. Feb 2019. doi: 10.1 10.10337hyySeeeec122200000002 URL https:://inkk aps org/doi/10 11t..p/eee.///yyeeeeeett.... 020602
Xiaoxia Wu, Ethan Dyer, and Behnam Neyshabur. When do curricula work? In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021. URL
httts://openreviiet net/ net/fore cortt//rrmmmae wamimmmiidttttt tW4QEEIppp
Ke Frédéric Boussemart, Fred Hemery, and Christophe Lecoutre. simple model generate hard satisfiable instances. In Leslie Pack Kaelbling and Alessandro Saffiotti, editors, IJCAI-05, Proceedings of the Nineteenth International Joint Conference on Artificial Intelligence, Edinburgh, Scotland, UK, July 30 August 5, 2005, pages 337-342. Professional Book Center, 2005, URL 0989. pdf]
Gal Ychuda, Moshe Gabel, and Assaf Schuster. It's not what machines can learn, it's what cannot teach. In Proceedings the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event, volume 119 Proceedings of Machine Learning Research, pages 10831-10841, PMLR, 2020. URL http::/pproceeddngg mlr. www.ppreessviive h html
Dinghuai Zhang, Hanjun Dai, Nikolay Malkin, Aaron C. Courville, Yoshua Bengio, and Ling Pan. the flows tell: Solving graph combinatorial optimization problems with gflownets. CoRR, abs/2305. 17010, 2023. doi: 10 200 7010. URL battps://dot org/10. 10. 48550/ 48550/arXiv arxiv. 2305. 17010
15

---

A Appendix
A.1 Evaluation Metrics
We use different metrics to evaluate the performance of models on CO problems. We assess the quality of an individual solution by the associated value of energy function (o) which represents the size of the solution sets in MIS and MVC. Here m refers to the problem instance under consideration. Optimal solutions Jopt are obtained with the Gurobi solver [Gurobi Optimization, LLC 2023]. For models that generate solutions indeterministically
sample ng different solutions oj per problem instance and calculate the test dataset average of the hest relative error Gel-
1 M
Em (app))p)) - Em
(9)
[mmp())))
where M denotes number of problem instances the test dataset. case deterministic algorithms only one sample is generated, i.c. ns = L. In several experiments it insightful to investigate the average relative error e,a for which we take the average instead of the minimum in Eq. 9 Analogously, we also define the best approximation ratio A A by
M Σ M
1
Em (ar,))
AR
min
(10)
Em co
The average approximation ratio AR is defined by taking the average instead of the minimum operation in Eq. [0] For VAG-CO calculate the average over the ordered greedy (OG) sampling method (See. App(A.7)
We always report the evaluation metric together with the standard error three different seeds except for the experiment RRG- 100 MIS, where only results on one seed are reported.
A.2 Linear Integer Program Formulations
Gurobi we formulate the CO Problems studied in this paper in Linear Integer Program Formulation (LIP)
(see Tab. 5
Problem Type
ILP formulation
ge (0. 13
MVC
min
Pi
2o + Øj if (1,5) E
MIS
min
qi
qi Φ i (1,j) 8
MaxCI
min
4:
9: 43 & 1 if (i.j) E
N Σ esj
max
MaxCut
ee(0.1)
(K.))EE
s.t.
Cij & q: + 9j if (i.j)
Gi < 2 (g. + 4;) if (1,))
Table 5: Integer Linear Program (ILP) formulations of MIS, MVC, MaxCl and MaxCut.
A.3 Ensuring Feasible Solutions
Since there is no rigorous guarantee that the model samples only feasible solutions that satisfy the constraints, we use fast post processing procedure to make sure that only feasible solutions are sampled. Here, we make use of our choice of the relative weighting of the energy terms A and B (see Tab, [T] in the Ising formulation, which ensures that only feasible solutions minima in the energy landscape. Therefore, can detect violations when Eg > and search for the spin causes the largest amount violations. Subsequently, we change the spin value of the node with the highest number of violations to satisfy the constraint and repeat the process until Eg O. We observe in our experiments that this post proccerocessing step is typically unnecessary. since only in rare cases violating solutions sampled.
16

---

A.4 Standardization of the Energy Scale
Since the energy scale of CO problems vary significantly, good choice of hyperparameters like the initial temperature To, learning and relative weighting between the policy and value loss vary between different CO problem instances. Therefore, we standardize the energy scale that makes the choice good hyperparameters easier. For this purpose we first express the binary energy function E(q) (see. Tab in terms of spins, i.c. as E(or) by substituting 9 with 2i+1 We then sample states for CO problem instances from the training set by using Random Greedy Algorithm (RGA, see App. [A.6.2). Since we use only few RGA steps (see App. (..62). this algorithm performs only slightly better than a completely random algorithm. From these states the mean energy ERGA and the standard deviation (ERG)) over the training dataset computed. Subsequently, we standardize the energy scale by:
B(O)= = E(o) ÉRGA
(1I)
std (
A.5
Baseline Methods
We compare VAG-CO with our own implementation of Mean Field Approximation (MFA) and Erdos Goes Neural (EGN), In both of these approaches we consider the case with and without anncaling. In both methods the probability of a state a factorizes into product of independent Bernoulli probabilities. Therefore, the state probability is given by:
Eo(a E)
11 pola. EE)
(12)
While EGN and MFA both the same probability distribution they differ in how the energy, or in of annealing the free energy, is calculated. This is described the following.
A.5.1 Mean Field Approximation
We estimate the energy by sampling states from the parameterized probability distribution and apply REIN- FORCE with variance reduction to update the network parameters. The network parameter gradients for one CO Problem Instance are then calculated with
AO(E)
E
[[E(r)
b) Vo lo E)]
where b = Ex--p((r(s)) [B(()) For a batch of CO Problem Instances A6(E) is averaged over multiple E.
A.5.2
MFA with Annealing
In MFA-Anneal we estimate the free energy instead and again update the network parameters with REINFORCE with variance reduction explained above.
A.5.3
EGN
Since consider cases, where the energy function is known in Ising formulation (see Sec. the expectation E. Ea-pala L [E(r] can written down in closed form by substituting a; Po (o, Karalias and Loukas 2020 Therefore, the loss function for one CO Problem Instance can be computed with
L(E) E'(p((o EE)
and the network parameters can be updated directly by computing the gradients of this loss function.
A.5.4 EGN with Annealing
In case EGN-Anncal [Sun closed form with the following formula:
2022 the entropy of the parameterized probability distribution is calculated
H(pe(a|E)) -- (x(() log + (1 pe(or(E)) log(1 pa(a; E)]
17

---

A.5.5 Conditional Expectation
To obtain good samples from the MFA approach in a deterministic way, we adopt the Conditional Expectation (CE) method [Raghavan] 1988] as descrihed in Karalias and Loukas (2020)
To introduce randomness into the CF-based solution generation procedure we initialize the random node features (which are processed by the GNN) of every node with random vector with six hinary entries that are drawn from uniform distribution. Unless stated otherwise, we follow the procedure in Wang and Li (023) and report the best CE result of eight different random node feature initializations.
A.6
Greedy Algorithms
A.6.1 Degree Based Greedy Algorithm
The Degree Based Greedy (DB-Greedy) algorithm [Wormald] 1995] is polynomial time algorithm for the MIS problem. The DB-Greedy algorithm works in the following way: At first all nodes are sorted according their degrees. Then, starting from the smallest degree the node chosen to be part the independent set. the next step this node, its neighhors and their corresponding edges are deleted from the graph. These steps are repeated until the graph is empty.
This algorithm can also be applied to the MVC problem by using the fact that the complement of an independent set is vertex cover. In other words, nodes in the independent set are excluded from the vertex cover and nodes that are not part of the independent set are included into the vertex cover.
A.6.2 Random Greedy Algorithm
The Random Greedy Algorithm (RGA) is a general approach that can be utilized for solving a broad range of CO problems that can mapped to the Ising model. Initially, the algorithm randomly samples spin values with uniform probability. Then, for a fixed number of iterations a spin is randomly selected and its value is changed if it decreases the energy value. We set the number of iterations to N - ngg where N is the number of nodes in the graph and na is the number of repetitions per node.
For the purpose of the standardization the energy scale (scc. App. A.4 we employ this algorithm with ny set one.
A.7 Study on Sampling Methods
In order to find out, how the best samples can be obtained with VAG-CO we study in Fig. [ three different sampling methods on the RRG- 100 MIS and RB 200 MVC dataset. Here evaluate e%o on the test dataset and plot over the number of samples n, that are used for each graph in dataset. We also add MFA and DB-Greedy results, where we show for MFA how the method improves when more solutions are sampled. Since DB-Greedy is deterministic algorithm, we draw a horizontal line that indicates the solution quality of the algorithm. To draw relation to results that are presented in Fig. Tgg[[(eet [Left, Middle) we also add vertical dotted line at nış- 8 samples. For VAG-CO denote the first method as sampling (S), where for each graph ns solutions are sampled according to the corresponding probability distribution. In the second VAG-CO sampling method called ordered sampling (OS). we use for each graph no different BFS node orderings and sample ns/noo states per graph ordering. Finally, we sample VAG-CO solutions with method called ordered greedy (OG), where we generate solutions greedily for each BFS ordering with no n.s. Results in Fig. show that the sampling strategy OG always outperforms all other sampling strategies that we proposed for VAG-CO. Additionally, we see that as we increase no in the OS sampling strategy, the performance improves consistently. Remarkably, DB-Greedy exhibits the best solution quality when MFA and VAG-CO are allowed only one sample (ng 1). However, DB-Greedy is outperformed by VAG-CO OG already with a modest amount of ns > 1 samples.
A.8 Experimental Details
In this section we provide additional details to experiments that are presented in Sec. 6
Parameter Checkpointing. With VAG-CO we checkpoint over the course of training the parameters that obtain the best * and the best Egel on the validation set. For testing we always use the checkpoint with the best éyel. except for the results in Fig. 3 when the sampling strategies (S) and (OS) used. Here use the checkpoints with the best elel-
Hyperparameter Tuning. For AFFe-An the learning rate and initial temperature are tuned on the validation
dataset of Enzymes MIS via a grid search. We considered learning rates (Ir E [5 x 10 10-4, 1 x 10 4 4.5 x 10 "D and initial temperatures (1% e [0.5, 0.25, 0.1, 0.05)). With that, the number of GNN layers (L e [6, 8, 12I) are tuned. Finally, the annealing duration (N (Vammeat) is increased until we observe that longer annealing does not lead
18

---

WAG 100 MIs Daleset
RB-200
22"4
D3 Greedy
WAS CO.
2X10-2
VAG MGA Arench MITL:
WAG-00:5 5
NBG-rreah: CE MPA CE
V&G CD:
VAG-CO DS 7
12000 20
90 50 200
155
75
Figure 3: Ablation different sampling strategies for VAG-CO on the RRG-100 MIS dataset (Left) and on the RB-200 MVC Dataset (Right). The averaged best relative error dee on the test dataset is plotted over the number of solutions ns per graph. Error bars indicate the standard error over the test dataset.
to improvements. For MFA, we tested on ENZYMES MIS the learning rates (Ir E [5 x 10 x 10-4,
10 50-5.1 10 10-9I and with that tuned the numher of GNN layers (T. E [6. 8, 10, 12]) and stop training when no observed. In MFA and MFA-Anneal, when training on the other
significant improvement on the validation datasets we started with optimal ENZYMES MIS parameters and further tested different learning rates (Ir), initial temperatures (Tb) and number of GNN layers (L) individually on each dataset. For experiments on the RRG-100 MIS and RB-200 MVC we make a more extensive hyperparameter search. In RRG-100 MIS for MFA- Anneal we first tuned the initial temperature in the range of To [0.25, 0.1, 0.05, 0.035, 0.025, 0.015, 0.01] and RB-200 MVC we search for the best initial temperature within the range of To [0.5, 0.25, 0.1, 0.05). For MFA and MFA-Anncal, we then iteratively search for the best learning rate in the range 10-1, 1,1x x 10 10-1, 5 x 10 10-1] and then for the best number of GNN layers L E [6, 8, 12]. In the experiments with EGN on RB-200 MVC we also performed extensive hyperparameter seach, where first the best temperature To was tuned within the range of To [0.25, 0.1, 0.05]. Then tuned the best learning rate within the range of (r [5- 10
E [5 x
10 1, 5. 10 -5, 10 .J Then we tuned for the best number GNN layers withing the range of L E [8, 10, 12, 14). In EGN without annealing the same hayperparameter search was performed but at Th O. The best hyperparameters of all methods are listed in A.1.
On VAG-CO we iteratively tuned the learning rate (Ir e |5 10-2, 10 initial temperature (1b [0.05, 0.7, 0.1]). number of GNN layers (T. E [3, 4]) and number of annealing steps (Nauneal E [3000, 6000) on the ENZYMES MIS validation dataset. We then initially test these hyperparameters on other datasets and adapt the number annealing steps and the initial temperature. The choice all hyperparameters in VAG-CO is listed in Tab.
Ablation on Subgraph Tokenization. In the ablation on subgraph tokenization in Fig. I we keep all hyperpa- rameters the same except for the time horizon T (see Sec. A.9 and the hyperparameter À in the PPO algorithm (see App. A.9.. For the subgraph tokenization run chose T
20, À 0.95 and l: 5. Therefore. the time horizon includes the generation of k 100 spins per graph during the data collection phase (see App. A.9 we would keep T the same, when we chose = 1 for the run without subgraph tokenization only T k 20 spins would be generated for cach graph. Therefore, we use T = 100 when no subgraph tokenization is used. As we always set X
4 this hyperparameter is adapted accordingly.
Experiments on Random Regular Graphs. Since GNNs suffer from node ambiguity on Random Regular Graphs (RRGs) Wang and Li 2023 using random node features Abboud et al. [20TT] can resolve this issue. Therefore, for the VAG-CO experiment in Fig. (L.ft) we sample standard Gaussian random node features with the dimension of six and concatenate them to the graph representation of VAG-CO. For MFA and MFA-Anneal random node features are already used (see A.5.5) and had not to be added for the RRG experiments.
Averaged results on hard instances. In Fig. [Left, Middle) only show results for specific generation parameters p on the RB-200 dataset and for specific values of d on the RRG-100 dataset. Here we report the corresponding averaged results on the RRG-100 MIS and RB-200 MVC dataset in Tab. VAG-CO significantly outperforms all other methods for ns 8.
A.9 Proximal Policy Optimization
Proximal Policy Optimization (PPO) [Schulman et al. 2017] is a popular RL algorithm that has two main components. The policy is represented by the network pa (o, G.) with parameters 0. The expected future reward is estimated by the value network Vo(G.) which is parameterized by d.
19

---

RB-200 MVC
RRG- 100 MIS
Evaluation metric A.1
AR*
(s/graph)
AR*
(s/graph)
DB-Greedy
1.0088 0.006
0.01
0.931
0.003
MFA:
1.0111
0.0002
0.10
0.917
0.06
MFA-Anneal: CE
1.0098 0.0001
0.10
0.955
0.06
EGN: CE
1.012 + 0.0001
0.10
N/A
N/A
EGN-Anneal: CE
1.0087 0.0002
0.10
N/A
N/A
VAG-CO (ours)
1.0059 0.0001
0.20
0.975
0,08
Gurobi (lmax 0.08)
1.0092 0.0001
0.06
0.945
0.07
Gurobi (fmax 0.1)
1.0073 0.0001
0.07
0.954
0.09
Gurobi (tmax 0.15)
1.0052 = 0.0001
0.09
0.974
0.12
Gurobi (lmax 0.2)
1.0037
0.0001
0.11
0.987
0.16
Gurobi (fmax 0.3)
1.0017
0.0001
0.13
0.995
0.22
Gurobi (fmax 0.5)
1.0006
0.0001
0.176
0.998
0.31
Gurobi (fmax = 1.)
1.0001 0.0001
0.26
1.000
0.37
Table 6: Best approximation ratio AR on the test dataset of the RB-200 MVC and RRG- 100 MIS dataset. Results where the AR* is closer to one is better.
Data Collection. In PPO the policy is rolled out through the environment in order to collect and store the states, actions, probabilities and the output of the value function into the rollout buffer. The data collection procedure completed after going through T steps, where T the so-called time horizon. Then, advantages
A, := A(o,, G.) and value targets V." ;= V (G.) are calculated by making use of the Generalised Advantage Estimation (GAE, Schulman et al. (0166). ), where A, is given by
A; = 6z (xA)6:+1 +... (72) (T+) 87-1.
(13)
with 6; = Ri + +V(++) - V(G.) and V A, V(+).. our experiments the reward R; given by
Eq. In our experiments we chose 7 = 1.0 and À = =1 +. We set T to he approximately equal to the average graph size in the dataset. During the data collection phase we always collect data from 1/ 30 different CO problem instances and for each instance sample nıs 30 solutions.
For the gradient updates in PPO, the loss is estimated with minibatches that randomly sampled from the rollout buffer. As the rollout buffer contains data size H ng we chose to sample data by first sampling H Haimi graphs and then for each graph, we chose to sample Nuunib solutions and for each solution Sminb time steps are sampled. We report the minibatch settings of cach experiment in Tab. 7
Rollout Buffer.
Training. As described before, during training minibatches of data Dauni is randomly sampled from the replay buffer. In PPO the overall loss that depends on the L(0, d; Dmintt) is given by L(O, &; Dminib) (0; Dminit) + ey La (á; Dminnt)). where the first term depends the policy and the second term on the value function. To specify a minibatch sample from Dminib, we will use the index nt.
The policy loss one minibatch sample is then defined as
(0).. = min (0) An, elippi An)
(14)
where the importance sampling ratio I, (0) paloniGai is used to compute a weighted Monte Carlo estimate
Pgr (rn)
of An. Here, 8' represents the old policy parameters that used fill the rollout buffer during the data collection phase and the new parameters that used in gradient descend. In Fq. [4 clipping function used to prevent that the probability of the current policy Po(a,IGr differs by more than a factor + from the old policy per (a.|GG.)
Similarly, the value function loss for sample is given by Le(d), (Vo(G.) (G.)))
The loss L(0, o; Dminit) is then updated for each minibatch in the rollout buffer so that each sample used least once. This procedure overall repeated Neepeat times, before we increment the epech step Napoch and new data is collected with the updated set of parameters. In our experiments we always chose = 0.1, cv 0.5 and Reepeat 2.
20

---

A.10 Model Details and GNN Architectures
Before processing the graph representation with GNNs, we encode its node features with an encoder MLP. To obtain node embeddings N, e E with E [L....:Z] LJ that depend on the graph structure, the encoded node features are processed by L layers of GNNs. Afterwards. graph embedkingg The global graph embedding is then concatenated to node embedding Nf where is the index of the spin whose value is to be generated (see Sec. (. We then use this embedding to calculate the policy and value function outputs (see App. A.9) with two separate MLPs.
apply global sum aggregation to compute a global
In our experiments, we employ two distinct message passing architectures, The first one is a Message- passiin Neural Network (MPNN) (Battaglia et al. 20188 but with lincar message functions. The second architecture also includes skip connections (MPNN-skip) Battaglaa et al. 20188 which use when more than three GNN layers are used.
Message Passing Neural Network. Our MPNN layer is defined by the following update for node embeddings:
N In
V Nr Σ Ea)) + Wakip NoC
(15)
JEN(()
where (ANG,,V,E, is a message update MI.P that computes messages with pairwise interactions between nodes, Skip connection are implemented by adding the term Watip after the node update, where Wakip E
Rd0+1)x4()) is weight matrix. Finally, layer normalization In[-] Ba et al.
2016 is applied it is commonly
done in Neural Network architectures when skip connections are used Liu el al. 2021].
VAG-CO we use a policy and value MLP with three layers. Both MLPs have 120 neurons in the first two layers. The value MLP has only one neuron in the output layer and the policy MLP has 2 output neurons, For the encoder MLP we use a two layer MLP with 40 neurons in cach layer. In case of the node update MLP y and the message update MLP & we 2 layers with 64 neurons in each layer. Layer normalization In| is applied after every layer within MLP except in the output layer of the policy and value MLP. We Rectified Linear Unit UU) activation functions Agarap (2018) except the output of the policy MIP,, where a softmax activation is used and except of the output of the value MI. where no activation is used.
In MFA with and without anncaling the model uses an encoder MI. with one layer and 64 neurons. The node update MLP v and the message update MLP & have the same number of neurals as in VAG-CO. The output MLP that is applied each node after GNN message passing steps the has three layers with 64 neurons each and the final output layer has neurons with softmax activation.
A.10.1 Subgraph Tokenization
For subgraph tokenization we have to make changes to the one-hot encoding as described in Sec. B and also to the model architecture as described in App./A.10) The one-hot encoding in the graph representation is adapted so that spins d that are going to he generated receive enumeration that indicates their position in the sliced list of BFS-ordered (sec Sec. 3) indices i + k. Then, the graph G, is processed by GNN that provides a node embedding GN No (ax) for each node. Along with a global sum aggregation, the node emheddings of the spins that are going to be generated are then concatenated according to their BFS order (Sec. (3) and further processed by the policy MLP that calculates the probability for each of the 2 spin configurations using a softmax output layer.
When the number of nodes N in the graph is not dividable by k. the number vertices in the CO problem instance description has to be increased without changing the inherent optimization objective. This can be realised by adding sufficient amount spins into the lsing model (see Sec. 2 with zero spin weight B and no connections to other spins.
A.11 Training Time and Computational Resources
All runs with VAG-CO were conducted either on A100 Nvidia GPU with 40 GB Memory or an A40 Nvidia GPUs with 48 GB Memory. In case of COLLAB MVC an A100 with 80 GB Memory is used.
The training time of our algorithm depends hy hyperparameters like the number of annealing steps N Nanneal (see 4.. the number of edges and nodes of graphs in the dataset, on the GPUs that used during training, on the time horizon and on the minibatch size that used for gradient updates PPO (see App. A.9) For
example the TWITTER MVC run with Nanneal = 4000, T = 30, Hmiuth = N Nuinih = Smint = 10 takes one
day. when trained on an 1100 40 GB Nvidia GPU. The run on ENZYMES MIS trained on a A40 GPU takes ten hours, when trained with N Nanesal 6000, T 20, Hmiii Nminbb 15 and Sminib 10.
21

---

A.12 Time measurement
Gurobi. All Gurobi results conductedon Intel Xcon Platinum 8168 @ 2.70GHZ CPU with 24 cores. All Problems expect for MaxCut formualted LIP and set the number of threads for Gurobi 26.
Learned Methods. For all learned methods (EGN, MFA, VAG-CO) present the time per graph if the algorithm are implemented in parallelized manner, where each sampling process runs on separate thread. always report runtimes by using the jit compiled performance of the Graph Neural Network.
A.13 Hyperparameters
VAG-CO Hyperparameters.
All hyperparameters that change across datasets are listed in Tab are specified in the corresponding section (c.g. App. A.
whereas hyperparameters that stay the same
A.14
Dataset
CO Problem Instance
learning rate
Nanneal
To
L
Naiini
H Hoinib
Swint
T
A
TWITTER
MVC
1 x 10
4000
0.05
3
10
10
30
5
COLLAB
MVC
1 10
12000
0.05
10
10
35
5
MIS
1 x 10
4000
0.05
3
10
10
10
30
5
MVC
1 10
2000
0.05
10
10
5
MIS
5 10-4
2000
0.05
10
10
10
5
MaxCI
5 10
6000
10
20
ENZYMES
MIS
5 10
6000
0.07
4
15
15
15
25
5
MaxCI
5 10
0.03
10
20
PROTEINS
MIS
5 10
10000
0.1
15
15
10
35
5
MUTAG
MIS
1 10
2000
0.05
10
10
10
7
RB-100
MVC
10
2000
10"2
5
10-3 0.0
10
10
20
5
RB-200
MVC
5 10
20000
0.05
3
15
10
10
30
RRG-100
MIS
5 x 10
20000
0.1
4
10
10
10
20
5
BA 200-300
MaxCut
1 x 10
7000
0.08
10
10
30
Table 7: Hyperparemeters that are used in VAG-CO on different datasets.
MFA Hyperparameters.
All runs with MFA used a hatch size H/ of 32 and we sampled as = 30 solutions. Furthermore, we used 8 GNN layers, and 6 random node features per node. In the RB-200 MVC experiment 10 GNN layers were used, For
used learning rate of 1 10 2000 annealing steps Nanneal and start temperature Th of 0.1 except for RRG-10 100 MIS and for RB-200 MVC. In RRG-100 MIS a To 0.015 is used and in RB-200 MVC To = 0.05. For MFA with and without annealing we trained for 2000 epochs and used learning rate of 5 x 10 except for ENZYMS MIS and RB-200 MVC, where learning rate of 1 x 10 is used.
MFA-Anneal
EGN Hyperparameters.
All runs with EGN on RB-200 MVC have batch size of H = 32. Here, in EGN-Anneal we use 12 GNN layers and in EGN without anncaling 10 GNN layers performes best. In EGN-Anneal the initial temperature of To 0.1 with the learning rate of 5 x 10- obtained the best results. With EGN without annealing the learning rate used. Both methods are trained for 2000 epochs.
A.14 Annealing Schedule
As described Sec. we change the temperature in the reward Eq. 6 according a predefined annealing schedule. During the warm-up phase Nwarmup epochs the temperature is held constant at the initial temperature To. Afterwards the following temperature schedule is applied:
Napoch Nwarmup
T(Nopoch) Tanneat (Nepoch)
cos 2 +
+1
(16)
Nanneal
Here, T Tanneall (Nepoch) ) is the gradually decreasing amplitude of the temperature oscillations:
To
Taaneal Nepoch)
=
N Napoch Naarmuup
(17)
+
Manncal
22

---

1.0
Temperature
0.6
0.4
0.2
200
620
800
1000
1200
Nepoch
Figure 4: Illustration of the cosine modulated annealing schedule. The plot depicts the annealing schedule with Nwarmup = 200 and Nauneal = 1000 steps. The vertical red line marks the end of the warmup phase.
Here is a scaling factor that determines the slope of T anneall Nepoch). The parameter À determines the numher temperature oscillations
schedule and Namneal the total number of epochs that follow phase. We use À 3, N Warmmup = 400 and e 6. The course of the annealing schedule is illustrated Fig. 4 One reason for the usage of cosine modulation function is rather practical, namely that T = is reached multiple earlier stage for given of
warmup
times during training, which allows assessment of the training success hyperparmeters. the absence cosine modulation, we found harder to assess final performance before the entire anncaling was finished. Similar periodic schodules for lcarning rates have been proposed as variants of Stochastic Gradient Descent [Lohchilovv and Hutter 2017
A.15 Derivations
A.15.1
Free-Energy Decomposition into Rewards
In the following we show that using the reward defined in Eq. 6 consistent with the goal of minimizing the
free-energy defined in Eq. 2
The right-hand side of Eq. contains the expectation of the energy E(or) and that is proportional to the entropy po po(or). For energy we obtain following decomposition into individual steps of solution
generation process (Sec. 3):
N Σ
N
N Σ
B(o)
Σ Joo,o. B,o,
Σ Juj, + B,
AE
(18)
By using autoregressive factorization the entropy also be decomposed in the following way:
((po(o|E)) - [os(olE) = - E
> log pa (a
(19)
E
log pa (asla E)
Therefore, we can use this decomposition by using the reward R.G., 3) = AF + 1 logp(ag F3) Eo-pp |,, (R,,,,,)),, maximizing this reward will then be equivalent to
By the relation F(0; 3, E) minimizing the free-energy.
A.15.2 Details on Remark
We first restate Corollary of [Dudik al. 2007 in our notation. In the context of this Corollary our energy function E can be regarded as a single feature and, therefore, we use their result for n = 1. In addition, we consider the case in which the samples S are drawn from the target distribution, i.e. æ = p(B*). Since the terms Boltzmann distribution and Gibbs distribution be used interchangeably obtain:
23

---

Corollary (Corollary of Dudík al. 2007D. Assume that E is bounded in [0. 11. Let 3 minimize Cp(3) + A|g with À = / Then with probability at least 1 ò for every Boltzmann distribution p(3).
(3) & (3) + |8| 2 (2(2)))
Now consider the result for the case / and subtract Lp(s*) (3*). Remark follows applying the definition the Kullback-Leibler divergence Dxz z to left-hand side. We note that an equivalent statement, however, with consideration the Rademacher complexity of E can be derived from Theorem 2 in Cortes el al. 2013]]
24


Reviews:
[{"review_id": "dCKeliGKot", "reviewer": "Reviewer_ojQm", "summary": "This submission proposes a novel unsupervised framework for solving graph combinatorial optimization problems. The proposed method is coined as VAG-CO, which autoregressively generates solutions to CO problems via annealing / reinforcement learning. Corresponding theoretical analysis is provided and numerical demonstration on various datasets are conduct. ", "questions": "See above.", "limitations": "See above.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- This paper is well written, and the technique is solid.\n- The proposed autoregressive generation is appropriate for solving graph CO problems whose solutions are a set over the nodes, which could be represented with a binary vector. \n- The authors use RL to optimize and carefully design an MDP to achieve that.\n- An open sourced implementation is also provided.\n- Solid technical analysis is provided.\n- Various experiments on both simulated and realistic benchmark and on both MIS and MVC tasks.", "weaknesses": "This submission is great. A few issues are listed as follows.\n- It seems the method is also applicable for other CO problems such as Minimum Dominate Set. Would it possible for the authors to elaborate about at least the possibility?\n- It would be also good to evaluate under some subset of the MIS benchmark (https://github.com/maxiboether/mis-benchmark-framework).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thoughtful and constructive review. In particular, we are glad that the reviewer appreciates among other things our theoretical analysis, that the paper is \"well written\", and most importantly, that our method is \"solid\". \n\n**Applicability to other CO problem types, like Minimum Dominate Set (MDS):**   \nOur method is designed to be directly applicable to problems that can be written in the form of an Ising Hamiltonian (Eq. 1). This class of problems is very broad and encompasses all of Karp’s famous 21 NP-complete problems (Lucas [2014]).   \nThe reviewer asked specifically about the MDS problem. This problem can indeed be formulated as an Ising-type CO problem. This is possible by using the reduction from MDS to the Set Cover problem (Kann [1992]). The latter can then be expressed as an Ising-type problem as shown in Section 5.1 in Lucas [2014].\n\n**Benchmarking on MIS-Benchmark:**   \nWe agree that using the MIS-Benchmark would be interesting. Nevertheless, we decided to focus on benchmarks from the most recent publications in the field, like Karalias et al. [2022] and Wang et al. [2023]. We added 3 new problem settings on 2 new CO problem types (Max-Cut and Max-Clique) works. For Max-Cut we significantly outperform the very recent work by Zhang et al. [2023] on the Barabási & Albert (BA) dataset with 200 - 300 nodes (BA 200 - 300). In case of the Max-Clique problem we benchmark on the ENZYMES and IMDB-Binary dataset and significantly improve upon the corresponding results in Karalias et al. [2022]. The results are shown in Tab. 4 and 5 in the rebuttal file.\n\n\nKann [1992], “On the approximability of NP-complete optimization problems”, Doctoral dissertation, Royal Institute of Technology  \nLucas [2014], “Ising formulations of many NP problems”, arXiv:1302.5843  \nKaralias et al. [2022], “Neural Set Function Extensions: Learning with Discrete Functions in High Dimensions”, arXiv:2208.04055  \nZhang et al. [2023], “Let the Flows Tell: Solving Graph Combinatorial Optimization Problems with GFlowNets”, arXiv:2305.17010  \nWang et al. [2023], “Unsupervised Learning for Combinatorial Optimization Needs Meta-Learning”, arXiv:2301.03116  "}, {"review_id": "ERPqSmFQE5", "reviewer": "Reviewer_fy9s", "summary": "This paper explores the use of autoregressive deep generative models for solving combinatorial problems.\n\nStarting from an optimization problem, the authors reformulate it using a Boltzmann distribution on the solution set.\nThis distribution can then be approximated with a parametric family, typically a neural network, following a variational approach (minimization of the free energy).\n\nWhile many previous works use the independent (or \"mean field\") approximation, it is suggested that building the solution iteratively is more expressive, and thus performs better.\nIn this case, the autoregressive model used is a graph neural network, which gives a conditional probability for the current spin given past choices.\nTo increase inference speed, two tricks are introduced: subgraph tokenization (making several decisions at once) and dynamic graph pruning (removing fixed vertices from the graph). \n\nTraining is achieved using a temperature annealing scheme which favors gradual concentration around global optima.\nThis annealing process enjoys an interpretation in terms of sample complexity and curriculum learning.\n\nNumerical experiments on the maximum independent set and the minimum vertex cover problem show the promise of this new method.", "questions": "L121: Why BFS in particular? Is there a way to make the order itself parametric?\n\nL160: Is there a way to perform subgraph tokenization that does not scale exponentially with $k$?\n\nL198: Are there other insights from curriculum learning that we could draw inspiration from? I'm not at all familiar with this literature\n\nL215: What if our problem provides no natural way to define partial rewards?\n\nL312: Why change the evaluation metric from approximation ratio to relative error? Not sure I grasp the difference", "limitations": "The authors do not discuss many limitations of their work, although they mention the additional complexity.\n\nSocietal impact is irrelevant here.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "### Originality\n\nThis paper describes a new combination of known ideas, coherently motivated by the specific goal of solving hard optimization problems:\n\n- formulation of optimization in the language of statistical physics\n- variational approximation of the Boltzmann distribution\n- iterative construction of a solution\n- RL framework with partial rewards\n- temperature annealing \n\nAs far as I can tell, the training tricks related to tokenization and pruning seem to be novel.\n\n### Quality\n\nThe proposed model makes theoretical sense, and rests on solid foundations of previous work.\nI particularly enjoyed the interpretation of annealing as curriculum learning: as the temperature gets lower, the learning gets harder.\n\nAs for the numerical experiments, they appear to be thorough and include a wide variety of algorithms.\nThe authors also took care to test on harder instances whenever nearly optimal solutions were too easy to reach.\nThey do however exhibit a few inconsistencies, which I point out below.\n\n### Clarity\n\nn.a.\n\n### Significance\n\nThis paper contributes to a wider discussion on the best variational ansatz for combinatorial optimization.\nComparing mean field with autoregression is an important first step, and hopefully the field can evolve towards more generic structured representations.\n\nSubgraph tokenization is also an interesting idea, which deserves to be pushed further.", "weaknesses": "### Originality\n\n### Quality\n\nA questionable algorithmic choice is the arbitrary order in which nodes are processed.\nIt isn't obvious why BFS ordering of the vertices makes more sense than other options.\nThe same goes for subgraph tokenization: why take the $k$ next vertices in the BFS order, even though they might be very far apart from each other?\n\nAs for the experiments, the realistic datasets of the first batch are seemingly too easy, since every method implemented by the authors (including the baseline DB-Greedy) shows near-perfect accuracy.\nThe results of Table 2 are very surprising to me, as both EGN and MFA are mean field methods with conditional expectation decoding, yet they exhibit wildly different performances on the minimum AR.\nOn the other hand, the average AR is only reported for new methods, not for the state of the art.\nI would welcome more detail from the authors on this point.\nIn any case, the benefits of autoregression as opposed to independence are not definitively proven by this series of benchmarks.\n\nThe synthetic datasets of the second batch are designed to be harder to solve, but the associated plots do not include other methods from the state of the art, which were analyzed on realistic datasets only.\nAgain, clarification from the authors would be very welcome.\n\n### Clarity\n\nThe notations are sometimes a bit hard to follow.\nSeveral algorithms included in the benchmark suite are only mentioned earlier in the text.\n\n### Significance\n\nDepending on the validity and fairness of the benchmarks, the results may be less significant than announced by the authors. \nI look forward to the rebuttal period to enlighten this aspect.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thoughtful and constructive review and would like to respond on various points mentioned in the weaknesses and questions.\n\n### Quality: \n\n**BFS:**  \nSee answer to question on L121 below.\n\n**“the realistic datasets of the first batch are seemingly too easy”:**  \nThese are datasets that are used in recent works in this field. We agree with the reviewer that these datasets might be too easy (see L304). For this reason we decided to benchmark also on synthetic datasets that are known to yield hard MIS and MVC problems.\n\n**“The results of Table 2 are very surprising to me, as both EGN and MFA are mean field methods with conditional expectation decoding, yet they exhibit wildly different performances on the minimum AR.”**  \nA possible reason might be the training method: in EGN the gradients are calculated directly from exact expectation values while our MFA estimates gradients via REINFORCE. The latter approach might be less prone to getting stuck in local minima due to an increased variance in the gradient estimation.  \nInterestingly, the recent work of Zhang et al. [2023] also states in App. D that EGN is worse than any of their baselines on their MIS benchmarks.\n\n**“On the other hand, the average AR is only reported for new methods, not for the state of the art.”**  \nThe purpose of the average AR metric in Table 2 & 3 is to compare the best achievable solution quality for VAG-CO with mean-field based approaches when they do not use conditional expectation. Hence, these results can be regarded as an ablation study in which no non-learned decoding method is used but where we simply report the average AR of the sampled solutions. The computational requirements of including all other under-performing (in terms of $AR^*$) methods in this ablation study would be sizeable. Consequently, we decided to consider only the best performing mean-field methods.\n\n**\"The synthetic datasets of the second batch are designed to be harder to solve, but the associated plots do not include other methods from the state of the art….\"**  \nSince these experiments are computationally expensive we conducted them only for the best performing methods in terms of $AR^*$ in Table 2 & 3. However, due to a request by reviewer rhwA we added the results on synthetic dataset for Gurobi with several time limits to the rebuttal file (Fig. 1).\n\n### Clarity:\n\n**“The notations are sometimes a bit hard to follow. Several algorithms included in the benchmark suite are only mentioned earlier in the text.”**   \nIf specific weaknesses of the notation are pointed out to us we would be happy to improve them.\nWe will introduce each algorithm in the experiment section of the updated manuscript.\n\n### Questions:\n\n**L121: Why BFS in particular?**    \nBFS will typically order the nodes such that the next $k$ spins are likely to be directly connected to already generated spins and will typically have several generated spins in their neighborhood. Consequently, the newly generated spins receive direct information on the already generated spins via message passing. For example, with depth-first search (DFS) one would expect that the newly generated spins have fewer already generated spins in the neighborhood. Consequently, in DFS there would be less information flow from the already generated spins to the spins that are to be generated next. This should lead to an increased probability of generating tokens that are sub-optimal once they become connected to the previously generated spins.  \nInvestigating this question empirically could be an interesting direction for future work.\n\n**L121: Is there a way to make the order itself parametric?**  \nOne can certainly make the order of the spins parametric, e.g. by sampling the order of spins from a probability distribution that is obtained via attention. This would be an interesting direction for future work.\n\n**L160: Is there a way to perform subgraph tokenization that does not scale exponentially with k?**  \nOne can reduce the number of possible tokens by masking out tokens that violate constraints of the CO problem (e.g. the independence condition in MIS). Whether this would result in a sub-exponential scaling depends on the specific CO problem instance.\n\n**L198: Are there other insights from curriculum learning that we could draw inspiration from?**  \nYes, one could increase the size of the graphs throughout training like in Lisicki et al. [2020].  \nMore generally, one could cast the task of generating suitable training graphs for curriculum learning in neural CO as a task for a separate agent like in “Teacher-Student Curriculum Learning” by Matiisen et al. [2017]. Here a teacher learns to generate the tasks (e.g. graphs in CO) for a student that learns to solve the corresponding CO problem.\n\n**L215: What if our problem provides no natural way to define partial rewards?**  \nWithout partial rewards the RL setting would be changed to an MDP with sparse rewards that are received once a complete solution was generated. This RL problem would be harder but our method could still be applicable.\n\n**L312: Why change the evaluation metric from approximation ratio to relative error?**   \nThe reason was to use in all three subplots of Fig. 1 metrics where lower is better. For the approx. ratio (AR) lower is better for MVC while for MIS higher is better. We agree with the reviewer that it might be better to not use the relative error in Fig. 1 (left) and to stick to $AR^*$. We will do so in the updated manuscript. The updated Fig. 1 can be found in the rebuttal file.\n\n\nZhang et al. [2023], “Let the Flows Tell: Solving Graph Combinatorial Optimization Problems with GFlowNets”, arXiv: 2305.17010  \nLisicki et al. [2020], “Evaluating Curriculum Learning Strategies in Neural Combinatorial Optimization”, arXiv: 2011.06188  \nMatiisen et al. [2017], “Teacher-Student Curriculum Learning”, arXiv: 1707.00183"}, {"review_id": "zKcHUM7KXm", "reviewer": "Reviewer_59Wh", "summary": "The authors present Variational Annealing on Graphs for Combinatorial Optimization (VAG-CO), a novel method for tackling combinatorial optimization problems. This method presumably combines elements of variational inference, annealing, and graph theory to form an effective optimization approach.\n\n=====\nAfter reading the rebuttal, the authors have convinced me to increase my score.", "questions": "How do the authors define hard problems in the paper? What are the criteria?\n ", "limitations": "The experiments are not extensive.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Novel Methodology: The paper presents a new approach to solve combinatorial optimization problems, Variational Annealing on Graphs for Combinatorial Optimization (VAG-CO), which is innovative and adds value to the existing literature.\n\nOvercoming Limitations of Existing Methods: The authors identify the limitations of the widely used Mean-Field Approximation (MFA) and propose a way to overcome these, which demonstrates a deep understanding of the problem space.\n\nImprovement in Efficiency: The use of sub-graph tokenization and entropy regularization to improve the efficiency of the training and inference process could potentially revolutionize the way these types of problems are solved.", "weaknesses": "1. The experiments are only completed on synthetic data. It will be more convincing if the proposed method can be applied and compared on some real data.\n2. Maximum Independent Set (MIS) and Minimum Vertex Cover (MVC) are well-known problems in graph theory and are widely used as benchmark problems in combinatorial optimization. However, they only represent a subset of the vast array of combinatorial optimization problems. A more extensive evaluation of the proposed method would involve its application to a larger and more diverse set of problems. This would enable a more comprehensive understanding of its capabilities and limitations. If the authors want to claim the wide application of the proposed method, they may want to conduct experiments for Traveling Salesman Problem, Knapsack Problem, Vehicle Routing Problem etc.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thoughtful and constructive review. We are happy to read that:\n- our method \"is innovative and adds value to the existing literature\",\n- we \"identify the limitations of the widely used Mean-Field Approximation (MFA) and propose a way to overcome these, which demonstrates a deep understanding of the problem space.\",\n- the \"The use of sub-graph tokenization and entropy regularization ... could potentially revolutionize the way these types of problems are solved.\".\n\nPlease find comments to the review in the following.\n\n### Weaknesses:\n\n**(1): “The experiments are only completed on synthetic data. It will be more convincing if the proposed method can be applied and compared on some real data.”**  \nWe are not sure about the reviewer’s definition of “synthetic”. Many of the graph datasets that we use like TWITTER, PROTEINS, COLLAB, IMDB-Binary, MUTAG, and ENZYMES are considered to be real-world datasets in the literature.\n\n**(2): “A more extensive evaluation of the proposed method would involve its application to a larger and more diverse set of problems.“**  \nWe managed to extend our experimental evaluation by adding benchmarks for MaxCut and MaxClique. For MaxCut we significantly outperform the recent work by Zhang et al. [2023] on Barabási & Albert (BA) graphs with 200 - 300 nodes. In case of the MaxClique problem we benchmark on the ENZYMES and IMDB-Binary dataset and significantly improve upon the corresponding results in Karalias et al. [2022]. These results are shown in Tab. 4 and 5 in the rebuttal file.\n\n### Questions:\n**(1): “How do the authors define hard problems in the paper? What are the criteria?”**  \nOur claim of hardness of MIS on Random Regular Graphs (RRGs) is based on the observation of a clustering transition for degrees $d > 16$ in the corresponding independent sets (see Barbier et al. [2013] for details). In the context of neural CO it was argued in Angelini et al. [2022] that these MIS instances are hard. They write:\n“We have argued that for the MIS using $d$-RRG with $d < 16$ is likely to be an easy problem and the test would be not very selective (we have in mind now only smart algorithms, not the GNN of Ref. [4] whose performances are so poor to be rejected even by easy tests). However, for larger $d$ we expect the optimization to become much more demanding because the clustering of the IS of large size is likely to create relevant barriers that affect any algorithm searching for the MIS.”. This claim is empirically compatible with our updated Fig. 1 left in the rebuttal file. It shows that all compared algorithms become worse as $d$ is increased from 3 to 20. Remarkably, if Gurobi is restricted to the same runtime as VAG-CO, i.e. 0.08 s (or even to a longer runtime of 0.1 s), it is outperformed by VAG-CO for $d > 3$.  \nThe hardness of the RB graphs for the MVC problem is based on a correspondence between the MVC and forced satisfiable SAT problems on RB graphs (see Xu [2004] for details). These graphs were also used recently in Wang et al. [2023]. In agreement with their observations we observe in Fig. 1 middle in the rebuttal file that for low values of the parameter $p$ all methods, including time-limited Gurobi, exhibit a rise in the approximation ratio, i.e. these problems become harder for all investigated methods.\n\nBarbier et al. [2013], “The hard-core model on random graphs revisited”, arXiv:1306.4121  \nAngelini et al. [2022], “Cracking nuts with a sledgehammer: when modern graph neural networks do worse than classical greedy algorithms”, arXiv:2206.13211  \nWang et al. [2023], “Unsupervised Learning for Combinatorial Optimization Needs Meta-Learning”, arXiv:2301.0311  \nXu [2004], \"BHOSLIB: Benchmarks with Hidden Optimum Solutions for Graph Problems\", http://vlsicad.eecs.umich.edu/BK/Slots/cache/www.nlsde.buaa.edu.cn/~kexu/benchmarks/graph-benchmarks.htm"}, {"review_id": "SlNtQzXZ4B", "reviewer": "Reviewer_rhwA", "summary": "This paper uses RL and Annealing to train a learned CO solver. They demonstrate that their algorithm empirically performs well over many problem types and dataset types. ", "questions": "1. During the RL process, how do you compute the reward of an intermediate graph where spins for certain nodes have not been assigned? Do you only compute the reward for nodes that have assigned spins?\n\nOverall, due to the limitations in the experiments, issues with experiments, and lack of contribution, I am voting to reject this paper. However, if the authors can run some further experiments comparing this paper with Gurobi and including the time taken by each algorithm, I believe this paper will be improved. I am flexible, so if the authors adequately address my concerns, I will increase my score accordingly. ", "limitations": "They have addressed the limitations.", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 1, "contribution": 2, "strengths": "1. This paper uses several important baselines.\n2. The breadth of the datasets chosen is good.", "weaknesses": "1. You are missing a comma in the first sentence.\n2. Use \\citep to put parentheses around your citations\n3. You misspell Reinforce in the table.\n4. The experiments seem relatively limited in that comparison is only for two problem types. Moreover, the authors do not compare with Gurobi to get a range on how solvers that are not learned perform on these problems. Also, the time taken by each algorithm is not used to compare these algorithms. Indeed, these metrics are important to see to compare.\n5. I do not understand the contribution of this paper. Indeed, it seems that this paper combines annealing from Sun et al. and RL for CO, which is also known. The subgraph solving seems to reflect many classical CO algorithms and is not new. The subgraph tokenization may be novel. In general, I'm a little confused about what the overall contribution of this paper is in the context of the literature.\n6. Also, the improvement in the solving ratio seems very small, enough such that hyperparameter tuning your algorithm and not tuning the baselines can make this difference.  \n\nI might also change the title of the paper. Currently, it closely resembles \"Annealed Training for Combinatorial Optimization on Graphs\" by Sun et. al.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thoughtful and constructive review and would like to respond on various points in the following.\n\n### Weaknesses:\n\n**(1-3) Latex citation style, comma, and typos:**  \nWe will update the manuscript accordingly. \n\n**(4)**\n  \n**Limited experimental evaluation on only two different problem types:**  \nWe managed to extend our evaluation by adding benchmarks for Max-Cut and Max-Clique. For Max-Cut we significantly outperform the recent work by Zhang et al. [2023] on the Barabási & Albert (BA) dataset with 200 - 300 nodes. In case of the Max-Clique problem we benchmark on the ENZYMES and IMDB-Binary dataset and significantly improve upon the results in Karalias et al. [2022]. The results are shown in Tab. 4 & 5 in the rebuttal file.\n\n**No comparison to Gurobi:**  \nAs pointed out in “A.1 Evaluation Metrics” all our metrics are relative to Gurobi solutions (references in L272 and L286 and in Table 2 & 3).\nWe will place a corresponding statement at the end of the first paragraph of the “Experiments” section. Importantly, our updated results in Figure 1 left and middle (see rebuttal file) show that VAG-CO outperforms Gurobi for similar time budgets on the hardest synthetic problem settings, i.e. for low values of $d$ and $p$.\n\n**Lack of runtime comparisons:**  \nWe agree that runtimes are an important aspect and added the corresponding times for all of our experiments. The rebuttal file contains runtimes for Table 2 & 3 and for all remaining experiments in Fig. 1. In Table 2 we currently cannot provide runtimes for the results that were reported from Karalias et al. [2022] since they did not provide this information for their MIS experiments.\nFigure 3 in the rebuttal file shows that Subgraph Tokenization does not only yield higher solution qualities (as already shown in Fig. 1, right) but also sizable runtime reductions.\n\n**(5) Relation to Sun et al. [2022] and RL for CO and unclear contribution and novelty:**  \nThe reviewer points out that Sun et al. [2022] also use annealing. This is correct, they use annealing in combination with a mean-field method. We cite them in the corresponding section on “Variational Annealing” in L249 (and in L36).  \nWe are also well aware that RL for CO is not at all new and we do not claim that is new. On the contrary, in the “Related Work” section in L211 we dedicate a subsection to “Unsupervised Learning and Reinforcement Learning” to discuss the relation of our work to RL for CO.\nWe would like to re-emphasize our main contributions, including:  \n\n(i) identification of a central limitation related to mean-field methods in numerous recent works,  \n(ii) introduction of Subgraph Tokenization to enable efficient autoregressive graph generation,  \n(iii) state-of-the-art results on several popular CO problems on real-world and on synthetic datasets,  \n(iv) a motivation for annealing from statistical learning theory.  \n  \nWith regard to the novelty the reviewer criticizes that “subgraph solving” is not new, which is correct and we do not  intend to claim that this is a novelty in any way. The reviewer acknowledges that Subgraph Tokenization “may be novel”. To the best of our knowledge it is novel. In case the reviewer still has significant remaining doubt we would like to understand the underlying reasons.\n\n**(6) Small improvements, hyperparameters could explain them:**  \nThe reviewer states that improvements with VAG-CO are rather small and that in such scenarios insufficiently tuning the hyperparameters of baselines could explain superior performance.  \nWe agree that on real-world graphs (Table 2 & 3) the achievable remaining improvements are indeed small. For this reason we add synthetic problems which are known to be hard. On these problems there is more margin with respect to optimal solutions. We compare MFA methods, VAG-CO, DB-Greedy. In the updated figures we added results for Gurobi with various runtimes and find that we can outperform Gurobi on the hardest problems. The improvements of VAG-CO with respect to the other methods on these problems are sizable and highly significant (Fig. 1, left and middle in the rebuttal file).  \nImportantly, we report our hyperparameter optimization strategies in detail in “A.7 Experimental Details” L630 ff. which shows that the hyperparameter search for the MFA methods in Fig. 1, left and middle was actually very extensive and certainly not less extensive than for VAG-CO.\n\nThe reviewer criticizes that the title is too similar to the title of Sun et al. [2022]. If it is possible we would follow the reviewer’s advice and change the title to e.g. “Beyond Mean-Field: Autoregressive Graph Generation for Combinatorial Optimization”.\n\n### Questions:\n\n**(1) How is the reward calculated for intermediate graphs?**  \nWe calculate the rewards only for the spins that have already assigned values (see Eq.6, L144, and  A.13.1 Free-Energy Decomposition into Rewards).  \n\nBased to the main points of this review we improved our work in terms of additional experiments, comparisons to Gurobi, and by adding runtimes. We hope that these improvements will be reflected in the updated score.\n\nZhang et al. [2023], “Let the Flows Tell: Solving Graph Combinatorial Optimization Problems with GFlowNets”, arXiv:2305.17010  \nSun et al. [2022], “Annealed Training for Combinatorial Optimization on Graphs”, arXiv:2207.11542  \nKaralias et al. [2022], “Neural Set Function Extensions: Learning with Discrete Functions in High Dimensions”, arXiv:2208.04055"}]

Author rebuttal:
We would like to thank the reviewers for their helpful reviews. Based on these reviews we could clarify several aspects of our work. In the following we will highlight the points of criticism that led to an extension of the presentation of our results and one point where we have difficulties in comprehending a fundamental point of criticism.

**Addition of runtimes**  
Reviewer **rhwA** rightfully argued that the work would be improved if we added runtimes for our experiments. We did so for all of our experiments and whenever this information was available for reported results. All runtimes are now provided in the updated figures and tables in the rebuttal file. The obtained runtimes show that our autoregressive VAG-CO exhibits similar runtimes as recent mean-field (MF) methods. The reason for this surprising result is that the MF methods rely on conditional expectation which is a time consuming procedure. Importantly, our autoregressive method would not be able to achieve such remarkable runtimes if we would not employ our newly introduced Subgraph Tokenization (ST) technique. As shown in Fig. 3 in the rebuttal file an increased $k$ for ST yields better results and substantially reduced runtimes. We thank reviewer **rhwA** for raising this question since it highlights the benefit of ST. 

**Experimental evaluation**  
We evaluate our method on 8 real-world datasets and two synthetic datasets for in total 14 different hardness settings. We report results for 13 different methods. Therefore, we are somewhat surprised that the extend of the experimental evaluation is criticized by reviewer **59Wh** (both weaknesses are related to this aspect) and reviewer **rhwA** (weakness 4: "The experiments seem relatively limited in that comparison is only for two problem types.").  
However, the reviewers are right in stating that our experiments focused on two CO problem types and that there are many others that would be interesting. Despite the very limited rebuttal period we, therefore, did our best to fully convince reviewer **59Wh** and reviewer **rhwA** of our experimental evaluation and can now report results on two additional CO problem types: Max-Cut and Max-Clique (Tab. 4 and 5 in the rebuttal file).  
For Max-Cut we demonstrate that VAG-CO outperforms the very recent work Zhang et al. [2023] on their version of the Barabási & Albert dataset with 200-300 nodes (BA 200- 300). Similarly, our new results for Max-Clique on the ENZYMES and IMDB-Binary datasets improve upon the results by Karalias et al. [2022]. In both cases VAG-CO represents the state-of-the-art. We include runtimes for all of these results and report the Gurobi performance with various runtimes.  
These additional results underpin the strong performance of VAG-CO and complement the numerous results that we already reported on MIS and MVC. In total, we are convinced that the updated experimental evaluation is extremely extensive in comparison to the standards in this field and that we could clearly establish the strong performance of VAG-CO.
  
In summary, VAG-CO is among the best performing methods on 10/11 real-world problems and the single best method on 8/10 real-world problems when no non-learned algorithmic components like conditional expectation are used. Since these real-world datasets that are frequently used in the recent literature are almost optimally solved we also include synthetic datasets that are known to be hard. VAG-CO outperforms the best learned methods on the real-world dataset on all 13 hard synthetic problem settings by a large margin.

**Unclear contribution and novelty**  
While three reviewers agree that this work represents a valuable contribution the remaining reviewer **rhwA** writes "I do not understand the contribution of this paper." and "In general, I'm a little confused about what the overall contribution of this paper is in the context of the literature.".  In the following we would like to discuss several aspects of this point and hope that we can resolve potential misunderstandings.  
Reviewer **rhwA** correctly pointed out several aspects of our work that are not new. As we argue below in our response to weakness (5) by reviewer **rhwA** this is correct but we never claimed these aspects to be new. For example, it is argued by the reviewer that the concept of annealing was already utilized in Sun et al. [2022] which  correct and also stated in L249. Then the reviewer calls into question the novelty of ST by writing that it “may be novel”. To the best of our knowledge it is novel and as our results show it is essential for the performance and runtime of VAG-CO. We hope that the reviewer will either outline the reasons for the doubt on this novelty or acknowledge it as an important new contribution.  
We would like to point out that our main contributions:

(i) identification of a central limitation related to MF methods in numerous recent works,  
(ii) introduction of Subgraph Tokenization to enable efficient autoregressive graph generation,  
(iii) state-of-the-art results on several popular CO problems on real-world and on synthetic datasets,  
(iv) a motivation for annealing from statistical learning theory,

are explicitly stated at multiple prominent places in the manuscript, including the abstract, the last paragraph of the introduction, and the conclusion. We hope that this clarification will facilitate the assessment of our contributions.  
  
Finally, we are convinced that all other points raised in the reviews were clarified in our rebuttals below and we are happy to address any questions in the discussion period.

Zhang et al. [2023], “Let the Flows Tell: Solving Graph Combinatorial Optimization Problems with GFlowNets”, arXiv:2305.17010  
Karalias et al. [2022], “Neural Set Function Extensions: Learning with Discrete Functions in High Dimensions”, arXiv:2208.04055  
Sun et al. [2022], “Annealed Training for Combinatorial Optimization on Graphs”, arXiv:2207.11542

Meta-review: The paper has found a majority of positive reviews, and nice formulation, some novel training tricks, significant but some conceptual limitations (processing order ...), solid technical method and analysis, good experiments were noted. On the downside, rhwA raised a number of concerns, which could be addressed to a significant degree. Overall, the AC follows the majority view of the paper.
The exhaustive rebuttal has been positively acknowledged.
Decision: Accept (poster)