Enhancing Column Generation by a Machine-Learning-Based Pricing Heuristic for Graph Coloring

Column Generation (CG) is an effective method for solving large-scale optimization problems. CG starts by solving a sub-problem with a subset of columns (i.e., variables) and gradually includes new columns that can improve the solution of the current subproblem. The new columns are generated as needed by repeatedly solving a pricing problem, which is often NP-hard and is a bottleneck of the CG approach. To tackle this, we propose a Machine-Learning-based Pricing Heuristic (MLPH)that can generate many high-quality columns efficiently. In each iteration of CG, our MLPH leverages an ML model to predict the optimal solution of the pricing problem, which is then used to guide a sampling method to efficiently generate multiple high-quality columns. Using the graph coloring problem, we empirically show that MLPH significantly enhancesCG as compared to six state-of-the-art methods, and the improvement in CG can lead to substantially better performance of the branch-and-price exact method.


Introduction
Branch-and-price is a widely-used exact method for solving combinatorial optimization problems (Barnhart et al. 1998) in the general form of Dantzig-Wolfe reformulation (Vanderbeck 2000). This formulation often provides a much stronger Linear-Programming relaxation (LP) bound than the more compact formulations of the same problem, which may lead to a significant reduction in the problem's search space. However, solving the LP can be challenging, because it typically has an exponential number of variables (or columns) that cannot be considered all at once.
Column Generation (CG) is an iterative method for solving large-scale LPs. As illustrated in Figure 1, CG starts by solving a subproblem with a small fraction of the columns in an LP, commonly referred to as the Restricted Master Problem (RMP). Then, the optimal dual solution of the RMP is used to set up a pricing problem to search for the column with the least reduced cost. If that column has a negative reduced cost, the column is included in the RMP to further improve its solution. Otherwise, the RMP has captured all the columns with non-zero values in the optimal solution of the original LP. Since an optimal LP solution typically has only Copyright © 2022, Association for the Advancement of Artificial Intelligence (www.aaai.org). All rights reserved. a small proportion of columns with non-zero values, CG is expected to solve the LP to optimality without the need to explicitly consider all the columns (Lübbecke 2010).
In the CG approach, repeatedly solving the pricing problem is typically a bottleneck (Lübbecke and Desrosiers 2005), because the pricing problem is often NP-hard. To attain computational advantage, heuristic methods are often preferred to trade column quality for computational efficiency. An exact method is usually used only if a heuristic method failed to generate any column with a negative reduced cost (Lübbecke 2010). Existing studies have explored a variety of pricing heuristics, such as greedy search (Mehrotra and Trick 1996;Mourgaya and Vanderbeck 2007) and metaheuristics (Taillard 1999;Malaguti, Monaci, and Toth 2011;Beheshti and Hejazi 2015). Based on past computational experience, Lübbecke (2010) notes that including multiple columns to the RMP at an iteration of CG can often speed up the progress of CG.
In this paper, we propose a novel Machine-Learningbased Pricing Heuristic (MLPH) for efficiently solving pricing problems. Specifically, we train an ML model offline using a set of solved pricing problems with known optimal solutions. For an unseen pricing problem at an iteration of CG, we use this ML model to predict the optimal solution of the pricing problem, which is then used to guide the search method to generate high-quality columns. To gain efficiency, we employ a linear Support Vector Machine (Boser, Guyon, and Vapnik 1992) for prediction and a sampling method for generating columns. As our method can potentially generate many columns, we introduce several column-selection strategies to form the new RMP to start the next iteration of CG.
By harnessing the knowledge learned from historical data, our MLPH has several advantages over ex-isting pricing methods: (1) compared to sampling-based methods (Dorigo, Birattari, and Stützle 2006;Cai and Lin 2016), MLPH can effectively generate columns with better reduced costs; (2) compared to other pricing methods (Gurobi Optimization 2018;Wang, Cai, and Yin 2016), MLPH can efficiently generate many more columns with negative reduced costs. As our MLPH can efficiently generate many high-quality columns, it can help CG capture the columns in an optimal LP solution with a fewer number of iterations.
We demonstrate the efficacy of our proposed MLPH on the graph coloring problem. Our experimental study shows that MLPH can significantly accelerate the progress of CG and substantially enhance the branch-and-price exact method.

Background
In this section, we first introduce different formulations of the Graph Coloring Problem (GCP). Then, we use GCP to illustrate the solving process of CG.

Graph Coloring Problem
GCP aims to assign a minimum number of colors to vertices in a graph, such that every pair of the adjacent vertices does not share the same color (Malaguti and Toth 2010). Let G(V, E) denote a graph, where V is the set of vertices and E is the set of edges. GCP can be formulated as: x i,c + x j,c ≤ z c , (i, j) ∈ E; c ∈ C, (3) x i,c ∈ {0, 1}, i ∈ V; c ∈ C, (4) z c ∈ {0, 1}, c ∈ C.
The binary variable z c denotes whether a color c ∈ C is used to color the graph vertices; and x i,c denotes whether a certain color c is used to color the vertex indexed at i. Since this formulation has a polynomial number of variables and constraints, it is commonly referred to as the compact formulation of the GCP, i.e., GCP-compact. Given that vertices with the same color must be part of an independent set, GCP-compact can be expressed as using a minimum number of Maximal Independent Sets (MISs) to cover all the vertices in a graph such that every vertex is covered at least once (Mehrotra and Trick 1996), which can be done systematically using Dantzig-Wolfe decomposition (Vanderbeck 2000;Vanderbeck and Savelsbergh 2006). The reformulated problem is commonly referred to as the Set Covering formulation of the GCP (GCP-SC), defined as: s.t. s∈S,i∈s x s ∈ {0, 1}, s ∈ S.
The binary variable x s indicates whether a MIS s is used to cover a graph, and S is the set of all the possible MISs in that graph. While GCP-SC provides a much stronger LP than GCP-compact (Mehrotra and Trick 1996), it can contain an exponential number of variables (or columns) to represent all the MISs in a graph. Hence, solving the LP of such a large-scale problem is challenging.

Column Generation
Given the LP of GCP-SC, CG aims to capture the columns with non-zero values in the optimal LP solution, starting from a RMP with a tiny fraction of the columns in the original LP: s.t. s∈S,i∈s Note that the integer constraints on x s are relaxed, and only a small number of MISs is considered initially, i.e., S ⊂ S. The RMP can be efficiently solved using the simplex method or the interior point method (Dantzig 2016), and its optimal dual solution π = [π 1 , · · · , π |V | ] associated to vertices (i.e., Constraint (10)) can be used to set up a pricing problem, to search for new MISs with the least reduced cost: The binary variable v i denotes whether the vertex i is a part of the solution, i.e., a MIS according to constraints (13) and (14). Note that the pricing problem for GCP-SC is the NPhard Maximum Weight Independent Set Problem (MWISP), where the weight of a vertex i is its dual solution π i to RMP. To tackle MWISP, related studies (Mehrotra and Trick 1996;Malaguti, Monaci, and Toth 2011) employ efficient heuristic methods. Only when a heuristic method fails to find any MIS with a Negative Reduced Cost (NRC), an exact method is used to solve the MWISP to optimality and so generate the MIS with the least reduced cost. If there exist NRC MISs, they are selectively included in the RMP to further improve its solution, according to a pricing scheme (Lübbecke and Desrosiers 2005). Otherwise, the RMP has captured all the columns in the optimal solution of the original LP, and hence the original LP is optimally solved.

Optimal Solution Prediction
We train an ML model to predict the optimal MIS of the MWISP by solving a binary classification task. In our training data, a training example (f , y) corresponds to a vertex in an optimally solved MWISP instance, where f denotes the feature vector that summarizes the property of the corresponding vertex and y holds a binary value of 1 (or 0) indicating whether that vertex is in the optimal MIS (or not).
We make use of several features that characterize a vertex of the MWISP, including vertex weight, vertex degree, and the upper bound of a vertex (defined by the sum of weights of that vertex and the vertices that are not adjacent to it). In addition, we adopt two statistical features (Sun, Li, and Ernst 2021) to further enhance the expressiveness of the feature representation for vertices. Given a sample of randomly generated MISs (s ∈ S), the first statistical feature measures the correlation between the presence of a vertex i and the objective values of the sample MISs, where s k i is a binary value, indicating whether the vertex i is a part of the k th sample; o k denotes the objective value of that sample; s i and o respectively denote, the frequency of the vertex i being in a sample and the mean objective value across all samples. A vertex with a high correlation score indicates that this vertex is likely to appear in the highquality MISs.
The second statistical measure uses the rank r of the sample MISs with respect to their objective values, A vertex with a high ranking score indicates that this vertex appears frequently in the high-quality MISs.
To gain computational efficiency, we adopt Support Vector Machine with linear kernel (linear-SVM) (Boser, Guyon, and Vapnik 1992) to best separate the positive examples (i.e., vertices in the optimal MIS) and negative ones (i.e., vertices not in the optimal MIS) in the training data. For a vertex i in an unseen MWISP instance, the prediction d i ∈ R of the trained linear-SVM is the distance of this vertex, in the feature space, from the optimal decision boundary. This indicates how confidently the linear-SVM classifies this vertex as either in the optimal solution or not according to the signed distance.

Generating Columns via Sampling
Based on the ML prediction, we can build a probabilistic model to sample multiple high-quality MISs. To generate a MIS, we start with a set containing a randomly selected vertex from the graph, and then iteratively add new vertices into the set until no new vertex can be added. We compute the probability of selecting a vertex via the ML pre- a logistic function to re-scale the prediction of a vertex into the range of [0, 1], and C denotes the set of candidate vertices not already adjacent to any vertex selected. The normalized value σ(d i ) can be interpreted as the 'likelihood' that vertex is in the optimal MIS. Note that we sample MISs starting from a random vertex, so as to increase the diversity of the generated MISs, which has an impact on the solving time of CG.

Columns Selection
MLPH can potentially generate a large number of NRC columns for unseen MWISPs, and adding all of these to the RMP at early iterations of CG can slow down the solving process of the RMP in the successive CG iterations. However, selectively adding NRC columns may increase the chances of missing out the optimal columns. Therefore, we empirically investigate several strategies to form the RMP at the next CG iteration: 1) add-all. Add all the newly generated NRC columns to the RMP. 2) add-partial. From the newly generated NRC columns, select a proportion of them to add to the RMP in increasing order of their reduced cost.
3) replace-existing. From all the columns, sequentially select columns for the next RMP in the increasing order of their reduced costs, while maintaining the diversity of the set of selected columns by skipping columns too similar to those already added. The algorithm is outlined in the Appendix.
At one end of the spectrum, the strategy (1) adds all the columns to the RMP, resulting in the fastest growth of the size of the RMP. On the other hand, the strategy (3) replaces some of the columns in the current RMP with newly generated columns, and it can maintain a fixed number of columns in the RMP. Due to this restriction, CG may require more iterations to capture all the optimal LP columns. The strategy (2) is somewhere in between these two extremes, resulting in relatively slow growth of the size of the RMP.

Experiment Settings
Graph benchmarks and problem instance generation. We use standard Graph Coloring Benchmarks 1 . Given a graph with n vertices, the goal is to find the columns with non-zero values in the optimal LP solution of GCP-SC, starting from a RMP initialized with 10n randomly generated columns. We note that reducing the number of samples can affect the performance of CG negatively. Among 136 benchmark graphs, we remove those whose initial RMPs already contain all the optimal columns and whose initial RMPs cannot be solved by an LP solver within a reasonable time. For the remaining 89 graphs, we label 81 of them as 'small' and 8 of them as 'large', according to the computational time for solving their initial RMPs. For a graph, we can generate multiple RMPs by seeding the initial set of random columns, and these RMPs can be viewed as individual problem instances because solving them can result in different optimal dual solutions and hence different subsequent MWISPs. For training, we generate 10 instances on 10 small graphs with random seed s = 1314. For testing, we generate 24 instances on each graph using random seeds s ∈ {1, 2, · · · , 24}, resulting in a total number of 1944 small instances and 192 large instances.
Data collection and training. For each training instance, we run CG using an exact, specialized solver TSM  to solve MWISPs to optimality. The MWISPs with optimal solutions are recorded every five CG iterations up to the 25 th iteration of CG. In the training data, the statistical features are computed from a set of n MISs, randomly sampled uniformly, with all features normalized instance-wise . The parameters for training SVM are set to the default values of (Chang and Lin 2011), except that the regularization term for misclassifying positive training examples is raised to the ratio between the negative training examples and the positive ones. For tuning the parameters in the logistic function, we employ Bayesian Optimization (BO) (Snoek, Larochelle, and Adams 2012;Nogueira 2014). Specifically, BO treats the MLPH as a black-box and attempts 300 runs of MLPH using different sets of parameters in the logistic function to minimize the reduced cost of the best-found solution of the MWISP at the first CG iteration.
Compared pricing methods.
• Gurobi, a state-of-the-art commercial Mixed-Integer-Programming (MIP) solver (Gurobi Optimization 2018). Such a MIP solver is used as the pricing method for GCP by related work (Malaguti, Monaci, and Toth 2011). By default, Gurobi aims to solve a MWISP to optimality, hence it spends most of its computational time on improving the duality bound. In addition to this default configuration, we include another setting, Gurobi-heur, that focuses on finding feasible solutions. This is done by setting the parameters 'PoolSearchMode' to 2, 'PoolSolutions' to 10 8 , and 'Heuristics' to 95%. • Ant Colony Optimization (ACO), an efficient meta-heuristic that has been investigated for many combinatorial optimization problems (Dorigo, Birattari, and Stützle 2006). ACO maintains a probabilistic distribution during the solving process and constructs solutions by sampling from that distribution. We adopt the ACO variant as described in (Xu, Ma, and Lei 2007). • Specialized methods. Since the optimal solution of a MWISP is the same as that of solving the Maximum Weight Clique Problem (MWCP) in its complementary graph, we also include three state-of-the-art MWCP solvers: 1) TSM , an exact solver based on the branch-and-bound framework with domainspecific knowledge for tightening the dual bounds; 2) LSCC (Wang, Cai, and Yin 2016), a heuristic method based on Local Search; 3) Fastwclq (Cai and Lin 2016), a heuristic method that constructs solutions in a greedy fashion with respect to a benefit-estimation function.
Computational budgets, evaluation criteria, and other specifications. Table 1 shows the computational budgets for solving small and large problem instances, respectively. In addition to an overall cutoff time for CG, we also set a  cutoff time for solving the MWISP at every CG iteration.
In particular, the exact methods are also subject to the time limit and are evaluated as heuristic methods. Moreover, we set for each of the pricing methods an individual termination condition. For LSCC based on Local Search (LS), we terminate LSCC if it cannot find better solutions for 50n LS iterations. This is because, empirically, LSCC can find highquality solutions efficiently, and providing excessive computational time can hardly improve these solutions. For MLPH, ACO, and Fastwclq, we set the number of constructed solutions to 50n. For TSM and Gurobi, they may terminate early when solving a MWISP instance to optimality. We use 'addpartial' as the default column-selection method with the column limit n. When a pricing method fails to find any NRC column, TSM is used to solve the MWISP to optimality and the optimal column is added to the RMP to start the next iteration.
We evaluate the performance of CG using a certain pricing heuristic based on two criteria, the computational time of CG for solved problem instances and the objective value of the RMP (the lower the better) for unsolved problem instances. The latter measures how close the solution of the RMP is to the optimal LP solution and so reflects the progress that CG has made. When reporting the results, we will address CG using a certain pricing method in short, e.g., CG-MLPH.
During CG the RMPs are solved by the default LP solver of Gurobi (Gurobi Optimization 2018). The experiment is conducted on a cluster with 8 nodes. Each node has 32 CPUs (AMD EPYC Processor, 2245 MHz) and 128 GB RAM. Our code is written in C/C++ and is available online 2 .

Results & Analysis
Results for CG using different pricing methods. Figure 2 shows the solving statistics for small graphs. The left sub-figure shows that CG-MLPH can solve many more problem instances than CG using other pricing methods (with a given computational budget). The right sub-figure shows that CG-MLPH can make more substantial progress than the comparison methods over all test instances. Noticeably, CG-Fastwclq and CG-ACO have comparable performances, and they are better than the remaining methods we consider. To better understand the CG using different pricing methods with respect to graph characteristics, we report in Figure 3 the best method on individual graphs. Overall, we observe that MLPH is particularly suitable for CG on relatively larger and/or denser graphs. In contrast, Fastwclp, ACO, and TSM are only competitive on small and sparse graphs. Table 2 shows the results on large graphs for CG. Note that the cutoff time is extended to 8, 000 seconds, and the   : 81 small graphs labeled by CG with the winning pricing method. MLPH is the best pricing method for 52 graphs, followed by TSM for 12 graphs, ACO for 8 graphs, Fastwclq for 7 graphs, Gurobi for 1 graph, and LCSS for 1 graph.  compared pricing methods are parallelized on 4 CPUs. It can be seen that CG-MLPH achieves the best performance on every individual graph, and it outperforms the other methods substantially in most cases. In contrast, the other methods are only competitive for certain graphs. Notably, CG-MLPH can optimally solve some problem instances, while other methods cannot solve any instance to optimality.
MLPH as a competitive pricing method. For the set of 8 large graphs, Table 3 compares the performances of dif-  ferent pricing methods for solving the pricing problems in the initial CG iteration (to ensure the results used for compared methods are from solving the same set of MWISPs).
Comparing MLPH with other pricing methods, MLPH can find many more NRC columns. Furthermore, the quality of the best-found column by MLPH is highly competitive. Subsequently, CG-MLPH achieves the best performance overall (Table 2). In addition, it can be noted that CG-MLPH solves all problem instances on the graph 'wap01a' to optimality using an average number of 36.5 iterations, while CG methods based on other pricing methods cannot optimally solve any of these, having at a minimum, an average number of 43.7 iterations. For other pricing methods, we can also observe that the good performance of CG with a certain pricing method is often accompanied by finding many high-quality columns by that pricing method, such as ACO for dense graphs (e.g., 'C4000.5') and Fastwclq for sparse graphs (e.g., 4-FullIns 5). Similar observations can be also made from the results on small graphs as shown in the Appendix.
The results indicate that finding a large number of highquality columns can accelerate the progress of CG. In particular, our proposed MLPH can find a large number of highquality NRC columns, thereby helping CG obtain much better LP objective values for unsolved problem instances or spending many fewer CG iterations for solved problem instances.
Efficiency and effectiveness trade-off in column selection. As shown in the left of Figure 4, adding all the NRC columns generated at every CG iteration (i.e., 'add-all') results in the faster convergence of the LP objective. On the other hand keeping a fixed number of columns in the RMP, by replacing the existing columns already in the RMP with newly generated NRC columns (i.e., 'replace-existing'), tends to slow down the progress of CG. This shows that adding more columns can increase the chance of capturing the optimal LP columns. When measuring the progress of CG in wall-clock time as shown on the right, we observe that 'add-all' and 'replace-existing' are comparable because 'addall' increases the computational burden for solving the fastgrowing RMP. Compared to these two methods, adding a proportion (n in our case) of the best NRC columns ('addpartial') better balances the trade-off between efficiency and  Table 3: Results for different pricing methods solving the MWISP at the initial CG iteration for large problem instances. The first statistic shows the number of columns with negative reduced costs that a pricing method can find. The second statistic shows the reduced cost of the best column that a pricing method can find. Both statistics are averaged over 24 problem instances generated using the same graph. effectiveness.

Branch-and-price with MLPH
In this part, we use CG-MLPH to enhance Branch-and-Price (B&P), an exact method that solves a GCP to optimality by recursively decomposing the original problem (root node) into subproblems (child nodes). During the solving process, CG is used at every node to compute their LP bounds (lower bounds), and a node can be safely pruned without further expansion if its lower bound is no better than the current best-found solution.

Setup
We use the B&P code from an open-source MIP solver, SCIP (Gamrath et al. 2020). For CG, in particular, an efficient greedy search is used as the pricing heuristic for tackling MWISPs. Only when it fails to generate any NRC column, an exact method called t-clique is used to find the optimal column to either certify the optimality of the LP or start the next CG iteration with the optimal column added to RMP. Once the LP at the current node is solved, the node is branched into child nodes, and the columns generated at this node are passed into the child nodes. Furthermore, the B&P implementation incorporates specialized techniques to GCP from previous studies (Mehrotra and Trick 1996;Malaguti, Monaci, and Toth 2011), and more details can be found in the Appendix. We refer to the default setup of B&P as B&P-def, and compare it with B&P-MLPH that replaces the greedy search in B&P-def with the MLPH for solving MWISPs. Although the greedy search has negligible computational cost, it is less effective as it can only construct a single column at a CG iteration. In contrast, our MLPH can sample many high-quality columns effectively. Empirically, we examine the sample size λ of MLPH in {10n, n, 0.1n}, and observe that no single sample size can fit all graph benchmarks. To best contrast our B&P-MLPH with B&P-def, we report the results when λ = 10n. From newly generated NRC columns, we add at most θ columns into the RMP in the increasing order of their reduced costs, θ = n for the root node and θ = 0.1n for child nodes. We observe that setting a smaller θ in column selection for child nodes can reduce the memory required for storing all the columns generated during the B&P process without sacrificing the performance. This is be- cause child nodes often have a sufficient number of quality columns inherited from their parents and their initial RMPs are already close to the optimum. For each method, a total number of 1584 seeded runs are performed to solve the GCPs on a set of 66 graphs in the Graph Coloring Benchmarks. The excluded graphs are either too easy (both methods can solve them within 10 seconds) or too hard to solve (both methods cannot solve the LP at the root node) under the cutoff time of 8000 seconds. When the LP at the root node is solved to optimality, we report optimality gap, defined as Gap = 100% × upper bound−global lower bound upper bound . The upper bound is the objective value of the best-found solution and the global lower bound is determined by the smallest lower bound amongst the remaining open tree nodes. Figure 5 shows the number of runs for which B&P-def and B&P-MLPH can obtain optimality gap within a certain threshold value. We can observe that 1) B&P-MLPH (red) solves GCP within a certain optimality gap in more runs than B&P-def (blue); 2) B&P-MLPH solves GCP to optimality (Gap = 0%) in 482 runs, better than 444 runs by B&P-def; 3) B&P-MLPH obtains optimality gap (i.e., the LP at the root node is solved) on GCPs in 1524 runs, whereas B&Pdef obtains optimality gap in 1325 runs.

Results
Next, we report the numerical results on 36 benchmark graphs where the performance of the two compared methods are significantly different (according to the student's t-test with a significance level of 0.05). The results are grouped into Tables 4-6 based on their comparative performances.    Table 6: Results for graphs solved by the two methods in some runs. Table 4 shows the results for graphs that can be optimally solved by both methods in all runs. Here, B&P-MLPH uses less solving time than B&P-def on 9 graphs, and the speedup is substantial on graphs such as 'school1' (38×) and 'flat300 26 0' (8×). In contrast, B&P-def performs slightly better than B&P-MLPH only on 3 graphs. Table 5 shows the results for hard graphs not solved by any method. B&P-MLPH can still solve the LP at the root node for most graphs and runs. However, B&P-def fails to solve the LP at the root node for many graphs, resulting in no optimality gap for those graphs. Table 6 shows the results for the remaining graphs. B&P-MLPH can solve GCP to optimality on more graphs and in more runs, compared to B&P-def. Overall, B&P-MLPH significantly outperforms B&P-def on 26 out of 36 graphs. Apart from these promising results, our studies also show that MLPH can be better integrated into B&P based on certain conditions. Firstly, if the MWISP can be solved by an exact method efficiently on a graph, then it is not necessary to use MLPH (or any other pricing heuristic). In particular, on the set of 10 graphs where B&P-MLPH does not outperform B&P-def, the performance of B&P without using any pricing heuristic is on par with B&P-def's. Secondly, if the improvement of the RMP becomes very slow, i.e., the tailingoff effect (Gilmore and Gomory 1961), an exact method can be used occasionally to solve MWISP to optimality even if MLPH can still find NRC columns. This is because MLPH is likely to keep finding NRC columns, which prevents the execution of the exact method from finding the optimal column and computing the Lagrangian lower bound. This reduces the chance of an early termination of CG. When applying this condition, we observe improved results of B&P-MLPH on 12 graphs. The detailed results can be found in the Appendix.

Related Work
Machine Learning for combinatorial optimization has received a lot of attention in recent years (Bengio, Lodi, and Prouvost 2021). Existing studies have applied ML in a variety of ways, such as learning variable selection methods Gasse et al. 2019;Liu et al. 2020;Furian et al. 2021) or node selection methods (He, III, and Eisner 2014;Furian et al. 2021) for exact branch-and-bound solvers; learning to select the best algorithm among its alternatives based on the problem characteristics (Liberto et al. 2016;Khalil et al. 2017b); learning to determine whether to perform problem reformulation (Kruber, Lübbecke, and Parmentier 2017;Bonami, Lodi, and Zarpellon 2018) or problem reduction (Sun, Li, and Ernst 2021; Ding et al. 2020); learning primal heuristics aiming to construct an optimal solution directly (Khalil et al. 2017a;Kool, van Hoof, and Welling 2019); and learning to select columns for column generation (Morabit, Desaulniers, and Lodi 2021).
To our knowledge, we are the first to solve GCP using B&P with a ML-based pricing heuristic. Our work is in line with the concurrent studies (Quesnel et al. 2022;Morabit, Desaulniers, and Lodi 2022), which tackle routing problems using B&P where ML is used to enhance a labeling algorithm for pricing. Our work is also related to recent studies in predicting optimal solutions for combinatorial optimization problems (Li, Chen, and Koltun 2018;Sun et al. 2021;Ding et al. 2020). These solution-predictionbased methods focus on effectively finding a single best solution for a single problem instance, whereas our pricing heuristic aims to efficiently generate many high-quality solutions by solving a series of pricing problems.

Conclusion & Future Work
This paper presents a Machine-Learning-based Pricing Heuristic (MLPH) for tackling NP-hard pricing problems repeatedly encountered in the process of Column Generation (CG). Specifically, we employ Support Vector Machine with linear kernel to fast predict 'the optimal solution' for an NP-hard pricing problem, which is then adopted by a sampling-based method to construct many high-quality columns efficiently.
On the graph coloring problem, we demonstrate the efficacy of MLPH for solving its pricing problem -the maximum weight independent set problem. We demonstrate that MLPH can generate many more high-quality columns efficiently than existing state-of-the-art exact and heuristic methods. As a result, MLPH can significantly reduce the CG's computational time and enhance the branch-and-price exact method.
In future work, we would like to extend our MLPH method to other combinatorial optimization problems such as vehicle routing problems. Our overarching aim is to develop a generic ML-based pricing heuristic to speed up CG and branch-and-price for solving the Dantzig-Wolfe reformulation of combinatorial optimization problems. This document is organized as follows. We present the column-selection algorithm, 'replace-existing'; additional details about training our MLPH; a description about the B&P implementation provided by SCIP (Gamrath et al. 2020); a discussion on how to adapt MLPH to other combinatorial optimization problems. Then, we present additional empirical results and analysis for CG and B&P. For some abbreviations (used here and also in our main paper), their full terms can be found in Table 1.

The Column Selection Method:
'replace-existing' The 'replace-existing' method can keep a fixed size of the RMP throughout the process of CG, by replacing some of the columns already in the current RMP with some newly generated columns at each iteration of CG. We describe the algorithmic procedure as follows. We initialize an empty array to store the MISs (i.e., columns) selected for the next RMP. Each entry in this array is associated with a vertex in the graph, and an entry can store up to k MISs containing that vertex, where k is the ratio between the number of MISs and the number of vertices in the graph. To ensure the solution of the next RMP is no worse than the current one, we firstly add into the array the MISs with non-zero values in the optimal solution of the current RMP. Then, we sort the rest of the MISs in ascending order of their reduced costs, and attempt to add them into the array one by one. A MIS can be added into an entry if satisfying two conditions: 1) the MIS contains the vertex associated with that entry and 2) the entry has less than k MISs. We note that the second condition is for maintaining the diversity of the columns in the next RMP, which is of importance (Lübbecke 2010). Lastly, we use the MISs in the array to form the RMP to start the next iteration.

More Details of Training MLPH
In this part, we describe for MLPH the training graphs, the training time, and the learned parameters. We use ten small graphs randomly selected from the graph coloring benchmarks to train our MLPH method. These  graphs are "3-FullIns 4", "queen12 12", "1-Insertions 6", "mug88 25", "DSJC125.5", "flat300 20 0", "flat300 26 0", "DSJC1000.9", "DSJC250.1", and "queen11 11". We run CG on each of the graph and collect the training data during the CG process. More specifically, we use an exact method TSM  to optimally solve the pricing problems, and collect solved MWISPs with optimal MISs every five iterations of CG to ensure the diversity of the training dataset. This results in a training dataset consisting of approximately 20, 000 training examples. Using the training data, linear-SVM can be trained within one second, and its parameters are reported in Table 2. Then, we optimize the Logistic model p(d i ; β 0 , β 1 ) = 1 1+e β 0 d i +β 1 using Bayesian Optimization (BO) (Nogueira 2014). Using the ten graphs (listed above), BO runs CG-MLPH 300 times with different sets of the parameters β 0 and β 1 in the Logistic model for solving the MWISP at the initial CG iteration, and identifies the set of parameters, β 0 = 9.7750 and β 1 = 12.5564, that yield best performance. The BO process can take several hours to complete.

The B&P Implementation
We use the B&P code from SCIP (Gamrath et al. 2020). In the process of CG, the B&P uses an efficient greedy search (Mehrotra and Trick 1996) as the pricing heuristic and an exact method called t-clique for solving the pricing problems. Besides that, this B&P implementation incorporates some specialized techniques for solving GCPs from previous studies (Mehrotra and Trick 1996;Malaguti, Monaci, and Toth 2011), such as a branching  method that operates on the original compact formulation of GCP, a local search method that finds the initial primal solution of GCP (i.e., an upper bound) to be used as initial LP columns for RMP, an early branching method that may terminate CG early and safely by comparing Lagrangian lower bound of the current RMP with its objective value. In addition, the B&P can also benefit from the functionalities provided by the SCIP solver, such as primal heuristics for finding better feasible solutions and column-management algorithms for efficiently re-optimizing RMPs during the CG process.
It should be noted that we have made the following changes to the default B&P: 1) We used the formula for calculating the Lagrangian lower bound in (Malaguti, Monaci, and Toth 2011) instead of its default one, in order to ensure the exactness of the B&P, and 2) we do not enforce limitations on the number of simplex iterations (SCIP uses simplex methods to solve LP) and the maximum number of columns generated during the CG process, following previous work.

Adapting MLPH to Other Problems
In order to adapt MLPH to another combinatorial optimization problem, designing problem-specific features and modifying the sampling method will be required. Take the Vehicle Routing Problems as an example, where the pricing problem is the Resource-Constrained Shortest Path Problem (RC-SPP). We can make use of the problem-specific features (e.g, the length and different types of resource consumption of an edge) and our statistical features, and train a linear-SVM to predict whether an edge is a part of the optimal path. In the testing phase, the prediction of linear-SVM can be used to guide the sampling of high-quality paths (i.e., columns). Figure 1 shows the set of 81 benchmark graphs labeled as 'small'. It can be seen that most graphs have less than 1000 vertices and their densities are around or less than 0.5. The rest of the graphs either are relatively dense or have a relatively large number of vertices. The statistics of individual graphs are shown in Table 3 and 4 for completeness. We generate 24 instances on each of the benchmark graphs using different random seeds, resulting in a total number of 1944 instances. Using these instances, we test CG with different heuristic-pricing methods under a cutoff time of 1800 seconds, and the numerical results are reported as follows. Table 5 compares the performance of CG using different pricing methods. Overall, our CG-MLPH consistently outperforms other methods by a large margin, in terms of all three metrics. With respect to the total number of solved instances, CG-MLPH can solve 1096 instances as compared to 964 instances solved by the second-best method, CG-Fastwclq. In terms of the average objective value of the RMP (indicating how close the current RMP solution is to the optimal LP solution), CG-MLPH achieves a value of 12.934 whereas the second-best method CG-Fastwclq obtains a value of 13.494. For solving time, CG-MLPH achieves an average of 125.6 seconds, much better than 199.2 seconds by the second-best method CG-TSM. Further, we report in Table 6 that CG-MLPH achieves top-k (k ∈ {1, 2, 3}) performances on individual graphs many more times than other methods. CG-ACO, CG-Gurobi, CG-TSM, and CG-Fastwclq are comparable with each other, but their performances are less competitive to CG-MLPH's. Table 7 shows, for 54 graphs whose problem instances can be solved by CG using different pricing methods at least once, the number of iterations spent (averaged over solved instances) by CG with each pricing method. As can be seen, CG-MLPH spends the fewest number of iterations for 30 graphs (results highlighted in bold), and it is also very competitive for many graphs in the remaining 24 graphs. In other words, MLPH helps CG better capture columns in an optimal LP solution using fewer CG iterations than other pricing methods typically. This is because MLPH can find many more high-quality columns than other pricing methods, as evidenced by the results in Table 8. Comparing to the samplingbased method ACO, the quality of the best-found solution by MLPH is much better. Comparing to more effective methods such as Gurobi or TSM, MLPH can find many more NRC columns than those methods, and note that MLPH is still highly competitive in terms of the quality of the best-found column. Noticeably, CG-ACO and CG-Gourbi-heur are also competitive in terms of the iteration number (Table 7), because they can generate multiple NRC columns (Table 8).

Extended Results and Analysis for CG
Besides those summary statistics, we also include the re-    sults for individual graphs for reference: Table 14 and Table 15 show the results for CG using different pricing methods; Table 16 and Table 17 show the results for a pricing method for solving the MWISP at the initial CG iteration.

Extended Results and Analysis for B&P
In this part, we firstly demonstrate that we can better integrate MLPH with B&P based on some conditions. Then, we report parameter studies for B&P-MLPH. Table 9 shows that, on the set of 10 graphs where B&P-MLPH does not outperform B&P-def, B&P without using any heuristic-pricing method (i.e, only using the exactpricing method) can perform on par to or even better than B&P-def. In other words, the exact method is very efficient for solving MWISPs for these graphs. In such cases, it is not necessary to use a heuristic-pricing method (e.g., greedy search or MLPH). In practice, one can detect such graphs by running an exact method (with a small cutoff time) for solving the MWISP at the initial iteration of CG. Whether to use a pricing-heuristic method can then be determined based on the performance of the exact method on this graph. Table 10 shows that, for some benchmark graphs, B&P-MLPH can benefit from forcing the exact-pricing method to optimally solve MWISPs occasionally in the latter CG iterations. When enforcing this rule to B&P-MLPH, the performance of B&P-MLPH is significantly improved for a set of 14 benchmark graphs. However, its performance is also significantly degraded for another set of 14 graphs. These results indicate that, for certain graphs, finding the optimal solution by an exact method and computing the Lagrangian lower bound can early terminate CG. However, the exact method should be given a certain cutoff time, because it can introduce computational overhead and degrade the perfor-   mance of B&P-MLPH for some graphs. Next, we report parameter studies for B&P-MLPH. Table 11 compares the default setting of B&P-MLPH (used in our main paper) with a variant that uses a constant columnlimit parameter θ = n for column selection. Note that our default B&P-MLPH sets column-limit parameter θ = n for the root node and θ = 0.1n for other nodes. It can be seen that the performance of B&P-MLPH with the default parameter setting is better overall. The reason is of two folds. Firstly, the child nodes often have a sufficient number of quality columns inherited from their parents, and hence keep adding a large number of columns can be less useful. Secondly, setting a small column limit parameter for child nodes can slow down the growth in the size of RMPs and reduce the computational overhead for an LP solver to solve RMPs. Table 12 and Table 13 present the results for varying the sampling size {10n, n, 0.1n} in B&P-MLPH, which shows that no single sample size of MLPH fits all benchmark graphs. These results indicate that B&P-MLPH can be possibly further improved by configuring B&P-MLPH according to the characteristics of the problem at hand. More specifically, the sample size can be set adaptively using information collected from the solving process, such as the speed of the CG convergence when using different samples sizes. Further, selecting a suitable sample size for B&P-MLPH can be viewed as an algorithm selection problem (Liberto et al. 2016;. We leave these possible solutions to future work as building an advanced parameter selection model requires significant effort and is out of the scope of the current study.   18.0 (7) 6.6 (24)

(24)
N  The number of CG iterations spent averaged over optimally solved test instances (out of 24) and the number of optimally solved instances (in brackets). If a method cannot solve any instance for certain graphs, then 'N/A' will be placed. Note that there are 54 graphs whose problem instances can be optimally solved at least once by CG. It can be seen that MLPH helps CG capture optimal LP columns with a fewer number of iterations on 30 of the graphs, and the performance of CG-MLPH is also competitive for the remaining graphs with a few exceptions.   Table 9: Comparing B&P using exact-pricing method directly (B&P-exact) with B&P-def that uses a greedy search on the set of 10 graphs where B&P-def is better than B&P-MLPH. For the three graphs on the top, B&P-exact outperforms B&P-def substantially, whereas B&P-def is better than B&P-exact for the graph 'r1000.5'. The two methods are comparable for the remaining graphs. The results indicate that the exact-pricing method is very efficient for solving MWISPs and hence, it is not necessary to use a heuristic-pricing method for these graphs.  Table 10: Comparing the default setting of B&P-MLPH (used in our main paper) with a variant that forces the exact method to run even if columns with negative reduced costs are found by MLPH, named B&P-MLPH-forced. The exact method runs occasionally under a frequency of 5 CG iterations when the objective value of the RMP improves less than 0.05 in successive 5 iterations. The results show that this variant is better than the default B&P-MLPH on 14 graphs (Group 2) and worse on a different set of 14 graphs (Group 1). The results indicate that, for graphs in Group 2, finding the optimal solution by an exact method and computing the Lagrangian lower bound can benefit B&P-MLPH from the early termination of CG. However, the exact method should be given a certain cutoff time, because it can take too much computational budget and degrade the performance of B&P-MLPH for graphs in Group 1.  Table 11: Comparing the default setting of B&P-MLPH (used in our main paper) with a variant that uses a uniform column-limit parameter θ = n for column selection, named B&P-MLPH-uni. Note that the default B&P-MLPH sets different column-limit parameters for the root node (θ = n) and other nodes (θ = n). The results show that the default B&P-MLPH is significantly better than the compared variant on 6 graphs (Group 1); the compared one is better on 3 graphs (Group 2); the performances of the two methods are comparable for other graphs not presented here.  Table 12: Comparing B&P-MLPH using sample size 10n (used in the main paper) and the one using sample size n, named B&P-MLPH-n. Group 1 shows the graphs where using sample size 10n is better and group 2 shows the graphs where using sample size n is better. The performances of the two methods are comparable to the rest of the benchmark graphs. The results show that no single sample size can fit all graph benchmarks.  Table 13: Comparing B&P-MLPH using sample size 10n (used in the main paper) and the one using sample size 0.1n, named B&P-MLPH-0.1n. Group 1 shows the graphs where using sample size 10n is better and group 2 shows the graphs where using sample size 0.1n is better.
The results show that no single sample size can fit all graph benchmarks.    Table 16: Statistics of solving the MWISP at the initial iteration of CG. The first category shows the number of found columns with negative reduced costs. The second category shows the most negative reduced cost of the best-found solution for MWISP. The third category shows the mean reduced cost among the columns with negative reduced costs. For a graph, the statistics are averaged over 24 problem instances generated using that graph.