text:
\section*{References}

\section{Introduction}

Many real-world problems can be formulated as Combinatorial Optimization Problems (COPs), which are in general NP-hard to solve. For some well-studied COPs such as the Travelling Salesman Problem~(TSP), there exist well-established methods~\citep{applegate2011traveling}. However, these methods typically incorporate specialized techniques designed by human experts, based on problem-specific mathematical tools which cannot be readily applied to other COPs. In recent years, leveraging Machine Learning (ML) for combinatorial optimization has attracted a lot of attention, because ML can automatically learn from data and may improve existing optimization algorithms in a principled way~\citep{lodi2017learning, bengio2020machine}. In particular, ML models have been shown to be effective in predicting optimal solutions for COPs~\citep{li2018combinatorial,lauri2019fine,sun2019using,ding2020accelerating,zhang2020nlocalsat,abbasi2020predicting}.

These approaches typically assume easy instances of the problem under consideration have already been solved to optimality in advance. This makes it possible to train a ML model to learn the mapping from the feature representation of a decision variable to its value in the optimal solution, supervised by solved problem instances with known optimal solutions as class labels. Taking TSP as an example, a ML model can be trained to classify whether an edge is part of the optimal tour or not, based on the features (e.g., length) of that edge. The trained ML model can then be used to predict the optimal solution values for decision variables in an unsolved TSP problem instance. The `predicted solution' (or the ML prediction) is likely to be infeasible, however it can be utilized in multiple ways to boost a search method, e.g., to prune the search space that is highly unlikely to contain high-quality solutions  \citep{lauri2019fine,sun2019using,ding2020accelerating}; to seed stochastic local search methods \citep{zhang2020nlocalsat}; or to guide tree search methods~\citep{li2018combinatorial, selsam2019guiding, shen21learning}.

To our knowledge the aforementioned approaches are \emph{single-shot}, i.e., using a ML model to predict the optimal solution for a COP only once, followed by repairing the ML prediction using a search method. Nevertheless, such \emph{Single-Shot} Solution Prediction (SSSP) approaches can fail to predict the optimal solution to sufficient accuracy, possibly resulting in failure in boosting the downstream search method. This is possibly due to the high complexity of NP-hard problems, particularly problems that are highly constrained and/or contain multiple optimal (or near-optimal) solutions.

To remedy this ML prediction issue, we propose a framework that enables a ML model to iteratively predict the optimal solution for a COP in \emph{multiple shots}, namely the Multi-Shot Solution Prediction (MSSP) framework. The high-level idea of our MSSP approach is illustrated in Figure~\ref{fig:mssp}. Given a trained ML model, MSSP uses it to predict the optimal solution of an unseen  COP instance to guide the search method, similar to SSSP. However, our MSSP can improve the prediction quality of ML models by harnessing the feedback from the search method, without re-training the ML model. Specifically, we make use of a set of statistical measures as features in a ML model~\citep{sun2019using}, which can extract useful information from the solutions found by the search method, and can inform the ML model which value for a decision variable is likely to be optimal. \added[id=yunzhuang]{This feedback mechanism makes an iterative improvement on the ML prediction possible. To see why, we highlight that the refinement of the ML prediction in the iterative process of MSSP comes from the update of statistical features, which are in turn based on the solutions found by the search method. As new and possibly better-quality solutions are found by the search method in each successive iteration of MSSP, more and more better-quality solutions can be used to approximate the solution space and provide more effective indications on the optimal solution.} In this study, we instantiate our MSSP framework with a Support Vector Machine~(SVM) \citep{boser1992training} as the ML model and a sampling method as the search method. However, this framework is more generic hence alternative ML and optimization methods can potentially be utilized.

\begin{figure}[t]
	\centering
	\begin{tikzpicture}[shorten >=1pt,auto,node distance=1cm, thick,main node/.style={color = marine, rectangle, fill=gray!10, rounded corners, inner sep=8pt}, scale = 1.0]
	\node[main node] (1) at (0,0) {ML Model};
	\node[main node] (2) at (6,0) {Search Method};

	\path[every node/.style={color = marine}]
	(1) [color = marine, line width = 0.30mm, ->, bend left = 15] edge node {predictions} (2);

	\path[every node/.style={color = red}]
	(2) [color = red, line width = 0.30mm, ->, bend left = 15] edge node {feedback} (1);

	\end{tikzpicture}
	\caption{Multi-shot solution prediction framework: the ML model's prediction of an optimal solution is used to guide a search method. The search method then provides feedback to refine the predictions of the ML model. This process is repeated iteratively.}
	\label{fig:mssp}
\end{figure}

\begin{figure}[t!]
    \centering
    \begin{tikzpicture}[baseline={(0,0)}]
    \begin{groupplot}[group style = {group size = 2 by 1, horizontal sep = 60pt},height=0.475\textwidth, width=0.475\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xmajorgrids=true, ymajorgrids=true,  log ticks with fixed point, major tick length=0.05cm, minor tick length=0.0cm, legend style={at={(0.18,-0.115)}, legend image post style={scale=1}, column sep = 1pt, legend columns=2,draw=none}, scatter/classes={%
	non-optimal={mark=*,blue}, optimal={mark=triangle*,red}}]
  \nextgroupplot[%
        title= \footnotesize Iteration \footnotesize 1,
    ]
	\addplot[scatter,only marks,%
    		scatter src = explicit symbolic,
		    visualization depends on=\thisrow{size}\as\size,
            scatter/@pre marker code/.append style={
            /tikz/mark size=\size}]%
    	table[x=x, y=y, meta=label, col sep=comma] {figures/feats_1.txt};
  \nextgroupplot[%
        title=  {\footnotesize Iteration \footnotesize 5},
        yticklabels={,,}
    ]
	\addplot[scatter,only marks,
    		scatter src = explicit symbolic,
		    visualization depends on=\thisrow{size}\as\size,
            scatter/@pre marker code/.append style={
            /tikz/mark size=\size}]%
    	table[x=x, y=y,meta=label,col sep=comma] {figures/feats_5.txt};
    \legend{ \footnotesize non-optimal, \footnotesize optimal}
    \end{groupplot}
    \end{tikzpicture}
    \caption{The benefit of employing \textit{multi-shot} solution prediction method using SVM on a TSP instance with 30 nodes. The axes are two statistical features described in Section \ref{subsection: features}, and each mark represents an edge in a TSP instance. The red marks denote the edges in the optimal (or shortest) tour, while the blue marks denote the non-optimal ones.}
    \label{fig:feats}
\end{figure}

Figure~\ref{fig:feats} illustrates the benefits of employing MSSP -- the decision variables can be better separated by performing more iterations of MSSP. Specifically, this visualization shows the localization of edges with respect to the two statistical features on a TSP instance with $30$ nodes. In the first iteration, it is difficult to separate the edges in the shortest tour (red marks) from those non-optimal edges (blue marks). In the subsequent iterations of MSSP, the statistical features are refined using the solutions found by the search method. After 5 iterations, the optimal edges have been shifted to the upper-right corner, while the non-optimal ones moved to the lower-left corner. As a result, the optimal edges are more distinguishable and easier to identify, and so is the optimal path.

To demonstrate the generality of our MSSP approach, we apply it to three well-known COPs, the traveling salesman problem, the orienteering problem, and the maximum weight clique problem. For each problem, we select for the SVM some problem-specific features and exploit the structure of the problem when sampling feasible solutions. We empirically show that, compared to existing SSSP approaches, our MSSP significantly improves the quality of the ML prediction and finds better-quality solutions. We further show that MSSP is very competitive compared to a range of other solution methods.

To demonstrate the wider applicability of MSSP, we use it as a pricing heuristic for column generation (CG)~\citep{lubbecke2005selected} and show that it can boost branch-and-price (B\&P)~\citep{barnhart1998branch} for solving the Graph Coloring Problem (GCP). Previous studies have attempted different heuristic methods for solving repeated NP-hard pricing problems in the computational process of CG for different COPs~\citep{lubbecke2005selected}. In our recent work~\citep{shen2022enhancing}, we have devised a SSSP approach and used it as a pricing heuristic to enhance CG and B\&P for solving GCP. Here, we extend \citep{shen2022enhancing} by using the newly developed MSSP as the pricing heuristic and show that MSSP can further improve the performances of CG and B\&P compared to existing pricing heuristics for GCP.

Our paper is organized as follows. Section~\ref{sec:rw} introduces recent developments in learning heuristics for combinatorial optimization. Section~\ref{sec:mssp-generic} presents the MSSP framework for tackling general COPs. Section~\ref{sec:mssp-specific} demonstrates how MSSP can be applied to three COPs, followed by computational results in Section~\ref{sec:ret}. Section~\ref{sec:cgbp} applies MSSP to enhance CG and the B\&P exact method. We conclude our paper in Section~\ref{sec:conclusion}.

\section{Related Work}
\label{sec:rw}

COPs can be expressed by minimizing a linear objective function of integer decision variables, subject to a set of linear inequalities~\citep{wolsey1999integer}. Since any integer variable can be represented as a string of binary variables, COPs can be formulated in the form:
\begin{equation}
\label{eq:cop}
P(\bm{v};\bm{c},\bm{A}, \bm{b}) \coloneqq \min_{\bm{v}} \{\bm{c}^T\bm{v}\;|\;\bm{A}\bm{v} \leq \bm{b},\;\bm{v}\in\{0,1\}^{n}\},
\end{equation}
where $\bm{c} \in \mathbb{R}^{n}$, $\bm{A} \in \mathbb{R}^{m\times n}$, and $\bm{b} \in \mathbb{R}^{m}$ are the problem data commonly referred to as cost coefficients, constraint coefficients, and constraint right-hand side, respectively.

COPs are typically NP-hard~\citep{daniel1994intro}. Heuristics can play an important role in finding good solutions efficiently~\citep{wang2016novel}. In the past decades, a variety of heuristics have been developed~\citep{zanakis1989heuristic,gendreau2005metaheuristics}, in particular metaheuristics, because of their versatility and ability in producing robust solutions to various COPs~\citep{sinclair1993comparison,blum2003metaheuristics,gendreau2005metaheuristics}.

Recently, learning heuristics via ML for COPs receives a lot of attention~\citep{bengio2020machine}. The existing methods can be classified into two classes: 1) one trains an agent via Reinforcement Learning (RL) and uses this RL agent to construct high-quality solutions sequentially~\citep{bello2016neural,nazari2018reinforcement,kool2018attention,deudon2018learning}; 2) the other trains a ML model under the supervision of optimal solutions and uses this ML model to predict the optimal solution that guides a search method to find high-quality solutions, i.e., SSSP approaches~\citep{li2018combinatorial,lauri2019fine,sun2019using,ding2020accelerating,zhang2020nlocalsat,abbasi2020predicting}. While RL-based methods are particularly effective for routing problems where a solution is an ordered sequence~\citep{kool2018attention}, SSSP approaches are widely applicable to general COPs~\citep{ding2020accelerating}. The latter is closely related to our work, as discussed in the following.

SSSP approaches collect the training data from a set of problem instances with known optimal solutions. Each training example corresponds to a decision variable in a problem instance, and consists of the optimal solution value of that decision variable as the class label and a feature vector containing its properties. Using these examples, a ML model can be trained to learn the mapping from the feature vector to the class label, i.e., to predict the optimal solution value of a decision variable based on its features. To obtain accurate predictions, existing studies attempt to derive a rich set of feature representations for decision variables. The most popular approach is learning feature representations using Graph Neural Networks~(GNNs)~\citep{li2018combinatorial, selsam2018learning, zhang2020nlocalsat, joshi2019efficient,ding2020accelerating,selsam2019guiding}. GNNs characterize a decision variable by incorporating the information from the variables in the local neighborhood, which is obtained naturally for graph-based problems~\citep{li2018combinatorial} or can be defined artificially according to constraint coefficient matrices for general COPs~\citep{ding2020accelerating}. Instead of learning feature representation using GNNs, \citet{sun2019using} devise statistical measures~(i.e., features) that directly characterize a decision variable in terms of its prevalence in a set of randomly sampled solutions and their objective values.

Given an unseen COP instance, a ML model trained for this COP can be used to predict the values of decision variables in the optimal solution. The ML prediction approximates the optimal solution and can be used to boost a search method in finding (near-)optimal solutions efficiently. Several approaches have been proposed to make use of a ML prediction, such as fixing a proportion of decision variables to their predicted values to reduce the problem size~\citep{sun2019using,ding2020accelerating}; rounding the ML prediction to integer values to seed the initial position of a stochastic local search method ~\citep{zhang2020nlocalsat}; using the ML prediction to guide a tree search method to firstly visit the nodes (i.e., sub-problems) likely to contain optimal solutions~\citep{selsam2019guiding,li2018combinatorial}; or using the ML prediction to build a probabilistic model for sampling~\citep{shen2022enhancing}. It can be observed that ML predictions have been successfully adopted to enhance generic methods~\citep{ding2020accelerating,sun2020boosting}, as well as specialized methods~\citep{sun2019using,zhang2020nlocalsat,selsam2018learning,li2018combinatorial}.

While existing approaches use a ML model to predict the optimal solution \emph{only once}, MSSP improves the prediction of a ML model by iteratively harnessing the feedback from a search method. This is done via the updating of the statistical features enabling the ML model to improve the prediction of the optimal solutions \textit{iteratively for  multiple times}. Noticeably, our proposed MSSP can be seen to be analogous to the popular deep neural networks~\citep{lecun2015deep} that update feature representations via a set of learned functions (or layers). At an iteration, the statistical features in the feature representations of decision variables are updated using the new solutions found by a search method, aiming to achieve a better ML prediction at the next iteration.  Moreover, MSSP can also be viewed as a generalization or variation of the estimation of distribution algorithms \citep{hauschild2011introduction}, as is discussed in Section~\ref{subsec:MHcompare} after introducing our method.
\section{Introduction}

Many real-world problems can be formulated as combinatorial optimization problems (COPs), which can be expressed by minimizing a linear objective function of integer decision variables, subject to a set of linear inequalities~\citep{wolsey1999integer}. COPs are generally NP-hard, hence solving them is challenging~\citep{daniel1994intro}. For some COPs such as the traveling salesman problem~(TSP), there exist well-established methods~\citep{applegate2011traveling}. However, these methods typically incorporate specialized techniques designed by human experts, based on problem-specific mathematical tools which cannot be readily applied to other COPs.

In recent years, leveraging machine learning (ML) for combinatorial optimization has attracted a lot of attention, because ML can automatically learn from data and may improve existing optimization algorithms in a principled way~\citep{lodi2017learning, bengio2020machine}. In particular, solution-prediction-based approaches have been developed~\citep{li2018combinatorial,lauri2019fine,sun2019using,ding2020accelerating,zhang2020nlocalsat,abbasi2020predicting}, where a ML model is used to predict the optimal solution for a COP, and the `predicted solution' (or the ML prediction) is then used to guide a search method to quickly find high-quality solutions.

For solution-prediction-based approaches, prediction of the optimal solution to sufficient accuracy is critical. However, this is a very challenging task to ML partially due to the high complexity of COPs~\citep{bengio2020machine}, particularly problems that are highly constrained and/or contain multiple optimal (or near-optimal) solutions. A common feature of these approaches is that they make only a single prediction of the optimal solution using their ML model, and such \emph{single-shot} solution prediction (SSSP) approaches can be inadequate in tackling hard COPs.

\begin{figure}[t!]
	\centering
	\begin{tikzpicture}[shorten >=1pt,auto,node distance=1cm, thick,main node/.style={color = marine, rectangle, fill=gray!10, rounded corners, inner sep=8pt}, scale = 1.0]
	\node[main node] (1) at (0,0) {ML model};
	\node[main node] (2) at (6,0) {Search method};

	\path[every node/.style={color = marine}]
	(1) [color = marine, line width = 0.30mm, ->, bend left = 15] edge node {Prediction} (2);

	\path[every node/.style={color = red}]
	(2) [color = red, line width = 0.30mm, ->, bend left = 15] edge node {Feedback} (1);

	\end{tikzpicture}
	\caption{Multi-shot solution prediction framework: the ML model's prediction of an optimal solution is used to guide the search method, which in turn provides feedback to refine the prediction of the ML model. This process is repeated iteratively.}
	\label{fig:mssp}
\end{figure}

This paper proposes a framework that enables a ML model to predict the optimal solution in multiple shots, namely \emph{multi-shot} solution prediction (MSSP).  The high-level idea of MSSP is illustrated in Figure~\ref{fig:mssp}. Given a COP instance, a ML model trained for this type of problem is used to predict the optimal solution to guide the search, similar to SSSP. However, our MSSP can improve the prediction quality of the ML model by harnessing the feedback from the search method.

\subsection{Literature Review}
Recent years have witnessed many studies leveraging ML to develop heuristic methods for COPs, which can be broadly classified into reinforcement-learning-based (RL-based) approaches and solution-prediction-based approaches. The former trains an agent via RL to construct high-quality solutions directly~\citep{bello2016neural,khalil2017learning,nazari2018reinforcement,kool2018attention,deudon2018learning,xinSCZ21}, which is particularly effective for routing problems with solutions being an ordered sequence. Compared to RL-based methods, solution-prediction-based methods are more applicable to general COPs because of their flexibility in the use of downstream optimization methods. The latter is closely related to our work, as discussed in the followi  ng.

Solution-prediction-based approaches typically assume easy instances of the problem under consideration have already been solved to optimality in advance. This makes it possible to train a ML model to learn the mapping from the features of a decision variable to its value in the optimal solution, supervised by solved problem instances with known optimal solutions as class labels. Taking TSP as an example, a ML model can be trained to classify whether an edge is part of the optimal tour or not, based on the features (e.g., length) of that edge. To obtain accurate predictions, existing studies attempt to enrich the feature representation for decision variables. The most popular approach is learning feature representation using graph neural networks~(GNNs)~\citep{li2018combinatorial, selsam2018learning, zhang2020nlocalsat, joshi2019efficient,ding2020accelerating,selsam2019guiding}. GNNs characterize a decision variable by incorporating the information from the variables in the local neighborhood, which is obtained naturally for graph-based problems~\citep{li2018combinatorial} or can be defined artificially according to constraint coefficient matrices for general COPs~\citep{ding2020accelerating}. Instead of learning feature representation using GNNs, \citet{sun2019using} devise statistical measures~(i.e., features) that directly characterize a decision variable in terms of its prevalence in a set of randomly sampled solutions and their objective values.

Given an unseen COP instance, a ML model trained for this problem can be used to predict the values of its decision variables in the optimal solution. The prediction of the optimal solution, if sufficiently accurate, can be used to guide a search method to quickly find quality solutions. Several examples include fixing a proportion of decision variables to their predicted values to reduce the problem size for an optimization method~\citep{sun2019using,ding2020accelerating}; rounding the ML prediction to integer values to seed the initial position of a stochastic local search method ~\citep{zhang2020nlocalsat}; using the ML prediction to guide a tree search method to visit first the nodes (i.e., sub-problems) likely to contain optimal solutions~\citep{selsam2019guiding,li2018combinatorial}; or using the ML prediction to build a probabilistic model for sampling~\citep{shen2022enhancing}. It is apparent that both generic methods~\citep{ding2020accelerating,shen21learning} as well as specialized methods~\citep{selsam2018learning,sun2019using} can benefit from solution prediction.

Besides developing heuristic methods, ML has also been used to accelerate the branch-and-bound (B\&B) exact method~\citep{land2010automatic} for deriving a certificate of optimality. In particular, recent studies~\citep{ quesnelWDS22,shen2022enhancing,abs-2201-02535} have applied solution-prediction-based approaches to boost column generation~(CG)~\citep{lubbecke2005selected} by solving a series of NP-hard pricing problems. The enhancement in CG can lead to substantial improvements in Branch-and-Price~(B\&P)~\citep{barnhart1998branch}, where CG is used to solve linear-programming relaxations during the computational process of B\&B. In addition, ML has been used to make critical decisions that may greatly influence the computational time of B\&B, such as approximating the strong branching decisions~\citep{lodi2017learning, gasseCFC019}, estimating the size of the B\&B tree to determine whether to perform a restart~\citep{hendel2021estimating}, deciding when to perform Dantzig-Wolfe decomposition~\citep{kruberLP17} or how to perform it~\citep{bassoCT20}, and others~\citep{vaclavikNSH18,morabitD021}. For a more comprehensive review, please refer to \citep{bengio2020machine}.

\subsection{Contributions}

This paper explores the idea of using feedback from search to further improve the prediction accuracy of a ML model. Our main contributions include the following:
\begin{itemize}
    \item We propose a novel framework, multi-shot solution prediction (MSSP). MSSP implements a feedback mechanism using statistical features, which extract useful information from solutions found by the search method and inform the ML model as to which value a decision variable is likely to be optimal. We design MSSP to maintain a set of best solutions sampled according to ML predictions across iterations. The evolving sample set can provide more reliable statistical information on the optimal solution, which can in turn help a ML model make more accurate predictions.

    \item We demonstrate the general applicability of our MSSP approach on three NP-hard COPs with different characteristics, including the maximum weight clique problem (MWCP), the orienteering problem (OP), and TSP. Empirical results show that MSSP improves the prediction quality of support vector machine substantially, resulting in competitive performance when compared with several classes of methods.

    \item We further demonstrate the efficacy of MSSP in boosting CG and B\&P with an application to the graph coloring problem, where MWCP is the pricing problem in CG and MSSP is used for solving MWCPs. Empirical results show that MSSP improves the performance of CG and B\&P compared to existing pricing methods.

\end{itemize}

\begin{figure}[th!]
    \centering
    \begin{tikzpicture}[baseline={(0,0)}]
    \begin{groupplot}[group style = {group size = 2 by 1, horizontal sep = 60pt},height=0.475\textwidth, width=0.475\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xmajorgrids=true, ymajorgrids=true,  log ticks with fixed point, major tick length=0.05cm, minor tick length=0.0cm, legend style={at={(0.05,-0.115)}, legend image post style={scale=1}, column sep = 1pt, legend columns=2,draw=none}, scatter/classes={%
	non-optimal={mark=*,blue}, optimal={mark=triangle*,red}}]
  \nextgroupplot[%
        title= \footnotesize Iteration \footnotesize 1,
    ]
	\addplot[scatter,only marks,%
    		scatter src = explicit symbolic,
		    visualization depends on=\thisrow{size}\as\size,
            scatter/@pre marker code/.append style={
            /tikz/mark size=\size}]%
    	table[x=x, y=y, meta=label, col sep=comma] {figures/feats_1.txt};
  \nextgroupplot[%
        title=  {\footnotesize Iteration \footnotesize 5},
        yticklabels={,,}
    ]
	\addplot[scatter,only marks,
    		scatter src = explicit symbolic,
		    visualization depends on=\thisrow{size}\as\size,
            scatter/@pre marker code/.append style={
            /tikz/mark size=\size}]%
    	table[x=x, y=y,meta=label,col sep=comma] {figures/feats_5.txt};
    \legend{ \footnotesize non-optimal, \footnotesize optimal}
    \end{groupplot}
    \end{tikzpicture}
    \caption{The benefit of employing \textit{multi-shot} solution prediction on a TSP instance with 30 nodes. The axes are two statistical features described in Section \ref{subsection: features}, and each mark represents an edge in the TSP instance. The red marks denote the edges in the optimal (or shortest) tour, while the blue marks denote the non-optimal ones.}
    \label{fig:feats}
\end{figure}

We provide an example to demonstrate MSSP and its benefits, using a TSP instance with $30$ cities. As illustrated in Figure~\ref{fig:feats}, the two sub-figures show the localization of edges with respect to the two statistical features at two iterations of MSSP. In the first iteration, it is difficult to separate the edges in the shortest tour (red marks) from those non-optimal edges (blue marks). In the subsequent iterations of MSSP, the statistical features of edges are refined based on the progressively improved solutions in the sample set. After 5 iterations, the optimal edges have been shifted to the upper-right corner, while the non-optimal ones moved to the lower-left corner. As a result, a ML model can better predict the optimal edges, and so is the optimal tour.

Our paper is organized as follows. Section~\ref{sec:mssp-generic} presents the MSSP framework for tackling general COPs. Section~\ref{sec:mssp-specific} demonstrates how MSSP can be applied to three tested COPs, followed by computational results in Section~\ref{sec:ret}. Section~\ref{sec:cgbp} applies MSSP to enhance CG and the B\&P exact method. We conclude our paper in Section~\ref{sec:conclusion}.

\section{Multi-shot Solution Prediction Framework}
\label{sec:mssp-generic}

The proposed MSSP framework aims to improve the prediction of a ML model by harnessing feedback from search iteratively. We first describe the key components of the framework, including the following: 1) a ML model that is trained to predict the optimal solution of a COP~(Section~\ref{sec:training}); 2) a sampling method that makes use of the prediction of the optimal solution to construct feasible solutions (Section~\ref{subsection: ps}); 3) a feedback mechanism from the search method to the ML model via statistical features~(Section~\ref{subsection: features}). Then, we outline the MSSP algorithm with the above components and discuss the design of MSSP in Section~\ref{subsec:mssp}.

\subsection{Model Training and Prediction}
\label{sec:training}
Following previous work \citep{li2018combinatorial,lauri2019fine,sun2019using,ding2020accelerating,zhang2020nlocalsat,abbasi2020predicting}, we construct the training data from a set of solved problem instances of a COP, each of which is associated with an optimal solution. Here, we consider only COPs with binary variables, but note that solution prediction can be also applied to general mixed-integer programs~\citep{abs-2012-13349}. In the training data, a training example corresponds to a decision variable $x_i\in \{0,1\}$ in a problem instance, which contains a feature vector $\bm{f}_i$ and a label $y_i$. The feature vector~$\bm{f}_i$ characterizes the properties of that variable, and the label $y_i$ holds its value in the optimal solution associated with that problem instance. We employ two types of features including the statistical features and the problem-specific features. We describe the statistical features in Section~\ref{subsection: features} and the problem-specific ones when introducing the test problems in Section~\ref{sec:mssp-specific}.

Using the training data, we can then train a ML model to best separate the training examples with labels $y=1$ and those with labels $y=0$, based on their feature representation $\bm{f}$. For the classification task, we employ support vector machine with the linear kernel~(linear-SVM)~\citep{boser1992training}, which aims to find the maximum margin hyperplane (parameterized by $\bm{w}$ and $b$) by solving the following quadratic programming problem with linear constraints:
\begin{align*}
	\min_{\bm{w},b, \bm{\xi}} & \quad \frac{1}{2}\bm{w}^T\bm{w} + r^+ \sum_{y_i=1}\xi^i + r^- \sum_{y_i=0}\xi^i, \\
	s.t. & \quad y_i\big(\bm{w}^T\bm{f}_i+b\big) \ge 1 - \xi^i, \quad i = 1, \cdots n_t,\\
	& \quad \xi^i \ge 0, \quad i = 1, \cdots n_t,
\end{align*}
where $r^+>0$ and $r^->0$ are respectively the penalty parameters for misclassifying positive ($y_i = 1$) and negative ($y_i = 0$) training examples; $\xi^i$ are the slack variables; and $n_t$ is the number of training examples.

After training, the linear-SVM model can then be used to predict the optimal solution values for decision variables in an unsolved problem instance. More specifically, for a decision variable $x_i$, we can compute the distance from its feature vector $\bm{f}_i$ to the optimal decision boundary of linear-SVM $(\bm{w}^*,b^*)$ in the feature space:
\begin{equation}
d_i=\frac{{\bm{w}^*}^{T}\bm{f}_i+b^*}{||\bm{w}^*||}.
\end{equation}
The signed distance $d_i$ indicates which value this decision variable is likely to take in the optimal solution. In an extreme case where $d_i$ approaches infinity, linear-SVM is confident to predict the optimal solution value of this decision variable to be~$1$. On the other hand, if $d_i$ approaches negative infinity, the optimal solution value of this variable is highly likely to be $0$. We then use a logistic function (parameterized by $\beta_0$ and $\beta_1$),
\begin{equation}
    p(d_i; \beta_0, \beta_1) = \frac{1}{1+e^{-(\beta_0 d_i + \beta_1)}},
\end{equation}
to normalize the distance value $d_i$ into a real number $p_i \in [0,1]$. For the unsolved instance, we can interpret the ML prediction $\bm{p} = [p_1, p_2, \cdots, p_n]$ of its optimal solution as the probabilities of the decision variables taking the value of $1$ in the optimal solution.

\subsection{Probabilistic Sampling as a Search Method}
\label{subsection: ps}

\begin{figure}[th!]
    \centering
\scalebox{1}{\begin{algorithm}[H]
\LinesNumbered
\SetAlgoLined
\KwIn{$\mathcal{P}(\bm{x};\bm{c},\bm{A}, \bm{b})$: problem data;\\
 \hspace{3em} $\bm{p}$: the ML prediction of the optimal solution.}

   $\mathcal{S} \gets \varnothing$, and $\mathcal{C} \gets \{v_1, \cdots, v_n\}$\;

   \While{$\mathcal{C}$ is not empty}{
        $i \sim \frac{p_i}{\sum_{j \in \mathcal{C}} p_j}$, where $i \in \mathcal{C}$\;
        $\mathcal{S} \gets \mathcal{S} \cup x_i$\;
        $\mathcal{C} \gets \mathcal{C} \setminus conflict(\mathcal{S})$\;
   }

    \Return{$\mathcal{S}$}
 \caption{Probabilistic Sampling.}
 \label{alg:ps}
\end{algorithm}}
\end{figure}

We make use of the ML prediction to form probabilistic models to sample solutions, namely the probabilistic sampling (PS) method. Our PS method is outlined in Algorithm~\ref{alg:ps}. Let $\mathcal{S}$ denote a partial solution containing decision variables fixed to $1$ and let $\mathcal{C}$ be comprised of unfixed variables (Line~1). A solution is constructed by iteratively fixing a variable $x_i$ to $1$, and the variable is sampled from a probability distribution over all the variables in the candidate set~$\mathcal{C}$ (Line~3). Here, the probability of sampling a variable~$x_i \in \mathcal{C}$ is proportional to its value~$p_i$ predicted by the linear-SVM, and the candidate variables are those that can be added to the current partial solution~$\mathcal{S}$. Once a new variable is fixed to~$1$, the conflict variables are detected and removed from the candidate set~(Line 5). Note that identifying and reducing the conflict domain of decision variables according to the already fixed decision variables is commonly referred to as domain propagation in the constraint-programming community, and generic techniques for domain propagation exist~\citep{AchterbergBKW08}. Finally, the process terminates and returns a feasible solution when all the variables are fixed, i.e., the set~$\mathcal{C}$ is empty (Line~2).

\subsection{Refining ML Prediction via Statistical Features}
\label{subsection: features}

We improve the prediction quality of a trained ML model (i.e., the linear-SVM model) by harnessing the feedback from search (i.e., the PS method). Specifically, we make use of two statistical measures~\citep{sun2019using} as a part of the feature representation for a decision variable. Given a sample of solutions from the PS method, these statistical features can characterize a decision variable in terms of its prevalence in the sampled solutions and their objective values. Once the statistical features for decision variables are updated, the ML prediction $\bm{p}$ can then be refined.

The first statistical feature computes a ranking score,

\begin{equation}
    f_r(x_i) = \sum^{K}_{k=1} \frac{\bm{s}^{k}_{i}}{r^{k}},
\end{equation}

\noindent where $\bm{s}^{k}_{i}$ denotes the value of the decision variable~$x_i$ in the $k^{th}$ solution in a sample set $\mathbb{S}$; $r^k$ denotes the rank of the $k^{th}$ solution in terms of its objective value; and $K$ denotes the sample size~$ K=|\mathbb{S}|$. A variable~$x_i$ with a high ranking score indicates that 1) the variable appears more frequently in the sampled solutions; and 2) the ranks of those solutions are high. Since this ranking-based feature accumulates the scores from sampled solutions, it is sensitive to the size~$K$. Thus, we normalize it using the maximum ranking score among all variables in a problem instance.

The second statistical feature computes the Pearson correlation coefficient between the values of a variable $x_i$ and the objective values of the sampled solutions,
\begin{equation}
\label{eq:corr}
    f_c(x_i) = \frac{\sum^{K}_{k=1} (\bm{s}^{k}_{i} - \overline{\bm{s}}_{i}) ({o}^{k} - \overline{o})}{\sum^{K}_{k=1} \sqrt{(\bm{s}^{k}_{i} - \overline{\bm{s}}_{i})^2} \sqrt{\sum^{K}_{k=1} ({o}^{k} - \overline{o})^2}},
\end{equation}

\noindent where $\overline{\bm{s}}_{i}$ denotes the mean of $x_i$'s values across the sampled solutions; ${o}^{k}$ and $\overline{o}$ are the objective value of $k^{th}$ solution and the mean objective value of the sampled solutions in $\mathbb{S}$, respectively. We normalize the correlation-based measure according to the type of the optimization problem. For a maximization problem, a variable highly positively correlated with the objective value is more likely to take the value of $1$ in high-quality solutions. Thus, this feature is normalized with respect to the maximum correlation score among variables. In a minimization problem, it is normalized by the minimum correlation score.

\subsection{Multi-shot Solution Prediction Algorithm}
\label{subsec:mssp}

\begin{figure}[th!]
    \centering

\scalebox{1}{\begin{algorithm}[H]
\LinesNumbered
\SetAlgoLined
\KwIn{$\mathcal{P}(\bm{x};\bm{c},\bm{A}, \bm{b})$: a problem instance; $\mathcal{ML}_{\bm{\theta}}$: a trained ML model;\\
      \hspace{3em}$T$: number of iterations; $M$: sample size.\\}

    $\mathbb{S}^{(1)}$: uniformly and randomly sampled solutions of $\mathcal{S}$\;
    \For{$t\gets1$ \KwTo $T$}{
      $\bm{F}^{(t)} \gets \bm{f}(\mathbb{S}^{(t)})$, where $\bm{F} = [\bm{f}^{(t)}_1, \cdots, \bm{f}^{(t)}_n]$\;
      $\bm{p} \gets \mathcal{ML}(\bm{F}^{(t)}; \bm{\theta})$, where
      $p_i = P(x_i = 1)$, $i=1..n$\;
      \For{$m\gets1$ \KwTo $M$}{
      $\mathcal{S} \gets$ PS($\mathcal{P}$, $\bm{p}$)  \tcp*{Alg. \ref{alg:ps}}

      Let $\hat{S}$ denote the worst solution in $(\mathbb{S}^{t})$\;
      \If{$\hat{\mathcal{S}} \not\in \mathbb{S}^{t}$ and $obj(\hat{\mathcal{S}}) < obj(\mathcal{S})$}{
            replace $\hat{S}$ with $S$ in $\mathbb{S}^{t}$\;
       }

      }
    }

    \Return{$\bm{p}$}
 \caption{\textit{Multi-shot} Solution Prediction Framework}
 \label{alg:mssp}
\end{algorithm}}
\end{figure}

Algorithm \ref{alg:mssp} outlines our MSSP method. At the beginning, a set of feasible solutions is generated uniformly and randomly (Line 1). Statistical features are then computed using the sampled solutions to form the feature representation of decision variables~$\bm{F}$~(Line 3). Based on $\bm{F}$, the trained linear-SVM with the logistic model is used to predict a `probability' value for each decision variable (Line 4). The ML prediction~$\bm{p}$ is then used to guide the PS method to construct $M$ solutions, and the unique and better-quality new samples replace the worst solutions in the sample set~$\mathbb{S}$~(Lines~5-9). The above process (Lines~3-9) is repeated for $T$ iterations, and the improved ML prediction is returned in the end (Line 10).

We design MSSP with an aim to \textit{construct and maintain a diverse set of high-quality solutions}. Firstly, maintaining a \textit{diverse} sample set allows MSSP to better explore the solution space and mitigate the possibility of premature convergence of the ML prediction to a local optimum. To this end, we employ the PS method to construct a diverse set of solutions, benefiting from the randomness in sampling. Secondly, continuous improvement of the \textit{quality} of solutions in the sample set is also vital, as the values of decision variables in high-quality solutions are likely to be highly correlated to the values in the optimal solution. Consequently, we replace worse solutions in the sample set with better-quality solutions. Finally, we note that MSSP may not always generate (near-)optimal solutions for certain COPs, as it aims to improve ML predictions steadily. More concretely, it has been shown that sampling-based methods (e.g., the PS method) are not sufficient to produce high-quality solutions for large TSP instances~\citep{joshi2019learning}. However, an improved ML prediction in MSSP can be used to boost an adequate optimization method to search for the best quality solutions, such as specialized search methods~\citep{selsam2019guiding,sun2019using}.

It is worth noting that MSSP is similar to the estimation-of-distribution algorithms (EDA)~\citep{hauschild2011introduction}, as they both alternate between randomly sampling the solution space and updating the probability distribution to increase the likelihood of producing `good' solutions. Their difference and connection are discussed in the following subsection.

\subsection{Relationship to Existing Metaheuristics}\label{subsec:MHcompare}

A fundamental difference that distinguishes MSSP from traditional EDAs is that, given a problem instance, EDAs fit a probabilistic model to a sample of solutions of that particular problem instance \textit{without supervision}, while MSSP leverages a ML model that is trained offline \textit{under the supervision} of optimal solutions for small problem instances. More concretely, Bayesian optimization algorithms~\citep{el2019updating}, a type of EDA, leverages Bayesian network to capture the relationship between decision variables. During the running of the algorithm, the Bayesian network is refitted at every iteration using the current set of sampled solutions. This refitting is very expensive and is one of the limitations of this type of algorithms. In contrast, the ML model in MSSP is trained offline using optimally solved small-scale problem instances and is not retrained during the problem-solving process.

MSSP can be naturally seen as a generalization of simple EDAs~\citep{dorigo1996ant}, i.e., those using a simple probability model and only updating the parameters rather than its structure. We demonstrate this connection with a well-known example, ant colony optimization (ACO)~\citep{dorigo1996ant}. This metaheuristic constructs a probabilistic model where the probability of setting a variable $x_i$ to one is proportional to $\tau_i^\alpha\eta_i^\beta$. $\tau$, the pheromone value, is essentially a statistical measure (similar to the statistical measures employed by MSSP) computed from sampled solutions; $\eta_i$ is a heuristic weight designed by human experts to guide the solution; $\alpha$ and $\beta$ are the parameters weighing these measures. While ACO relies on the two measures and combines them with a manually-defined function to form probabilistic models, our MSSP can incorporate many more useful measures (or features) and learn a model automatically via ML, under the supervision of solved problem instances with known optimal solutions.
\section{Multi-shot Solution Prediction for Three NP-hard COPs}\label{sec:mssp-specific}

In this section, we demonstrate how to apply MSSP to TSP, OP, and MWCP, respectively. For each COP, we present the problem formulation and describe the two problem-specific components of MSSP: the features in the linear-SVM and the algorithm for constructing feasible solutions in the PS method.

\subsection{Travelling Salesman Problem}
Given a set of $n$ cities, TSP aims to find a tour visiting each city exactly once such that the distance of the tour is minimized. Let the decision variable $x_{i,j}\in \{0,1\}$ denote whether the edge $(i,j)$ is in a feasible tour; let $d_{i,j}$ denote the distance of the edge $(i,j)$; let $u_{i}$ denote the visiting order of city $i$. TSP can be modeled using the Miller-Tucker-Zemlin formulation:
 \begin{align}
    \min_{\bm{x}} \;& \sum_{i=1}^{n}\sum_{j=1}^{n}d_{i,j}x_{i,j}, \\
        s.t. \;& \sum_{i=1}^{n} x_{i,k} = \sum_{j=1}^{n} x_{k,j} = 1;\; && 1\leq k \leq n,  \\
              & u_i - u_j + nx_{i,j} \leq n-1;\; && 2 \leq i,j \leq n,  \\
              & 0 \leq u_i \leq n-1;\; && 2 \leq i \leq n, \\
              & x_{i,j} \in \{0,1\};\; && 1 \leq i,j \leq n,

\end{align}
\noindent where Constraints~(6) specify that each city must be visited once and Constraints~(7) prevent sub-tours.

To apply MSSP to TSP, we train linear-SVM to predict whether an edge is in the shortest tour, supervised by a set of solved small TSP instances with known shortest tours (i.e., optimal solutions). Recall that the features for a decision variable (i.e., an edge) consist of statistical features and problem-specific ones. To initialize the statistical features, we use feasible tours randomly generated by permuting the order of cities. We derive problem-specific features to TSP from the length $d_{i,j}$ of an edge $(i,j)$, normalized by the mean or minimum length of edges connected to the city $i$ or $j$. Let $\overline{d}(i) = \frac{1}{n} \times \sum_{k=1}^{n} d_{ik}$ be the average length of edges connected to $i$.
Let $d^{max}(i)=\max_{k=1,\cdots,n} d_{i,k}$ and $d^{min}(i)=\min_{k=1,\cdots,n} d_{i,k}$ be the edge of maximum length and that of minimum length connected to $i$, respectively. The four normalized features are computed by
\begin{equation}
\frac{d_{i,j} - \overline{d}(i)}{d^{max}(i)-d^{min}(i)}, \;\;\;\;\;
\frac{d_{i,j} - d^{min}(i)}{d^{max}(i)-d^{min}(i)},
\;\;\;\;\;
\frac{d_{i,j} - \overline{d}(j)}{d^{max}(j)-d^{min}(j)}, \;\;\;\;\;
\frac{d_{i,j} - d^{min}(j)}{d^{max}(j)-d^{min}(j)}.
\end{equation}

Given an unseen TSP instance, the trained linear-SVM can predict whether an edge is likely to be in the shortest tour, and the ML prediction is used to guide the PS method for constructing new tours. Specifically, PS starts with the city indexed by $i \mod n$ ($i$ being the index to the current sample) and iteratively samples the next city to visit, given the current visiting city $i$ and past visited cities. The likelihood of visiting the city $j \in \mathcal{C}$ is proportional to the ML prediction $\bm{p}_j$ of the edge $(i, j)$. After all the cities are visited, PS returns to the first-visited city to form a feasible TSP tour. When updating the sample set~$\mathbb{S}$, we keep half of the initial random samples at the subsequent iterations. This maintains the diversity of the sample set~$\mathbb{S}$ and results in a more steady improvement of the ML prediction.

Since PS is not sufficient for the solving of large TSP instances~\citep{joshi2019learning}, we employ a greedy method to make use of the ML prediction for constructing high-quality tours, in addition to PS. At each iteration of MSSP, the greedy method constructs tours starting from different cities in a similar way to PS, except that the next visiting city (given the visited ones) is selected greedily rather than sampled randomly according to the ML prediction.

\subsection{Orienteering Problem}

OP aims to selectively visit a set of locations under some distance (or time) budget through a closed tour passing the depot location, such that the total prize collected from the visited locations is maximized. Let $i$ be the index to a set of $n$ locations associated with prizes~$p_i$; let $T$ denote the budget limit; let $d_{i,j}$ denote the distance of the edge~$(i,j)$; let $x_{i,j} \in \{0,1\}$ denote a decision variable that indicates whether the edge~$(i,j)$ is a part of the tour; let $u_i$ denote the visiting order of a location $i$.  Given the depot location indexed by $1$, OP can be formulated as
 \begin{align}
    \max_{\bm{x}, \bm{u}} \;& \sum_{i=2}^{n}\sum_{j=2}^{n}p_{j}x_{i,j}, \\
        s.t. \;& \sum_{j=1}^{n} x_{1,j}= \sum_{i=1}^{n} x_{i,1}=1, \\
              & \sum_{j=2}^{n} x_{k,j} = \sum_{i=2}^{n} x_{i,k} \leq 1;\; && 2 \leq k \leq n,  \\
              & u_i - u_j + 1 \leq (n-1)(1-x_{i,j});\; && 2 \leq i,j \leq n, \\
              & \sum_{i=2}^{n}\sum_{j=2}^{n} d_{i,j}x_{i,j} \leq T,\\
              & u_i \geq 0;\; && 2 \leq i \leq n, \\
              & x_{i,j} \in \{0,1\};\; && 1 \leq i,j \leq n,

\end{align}
\noindent where Constraint~(12) specifies that a feasible tour must start from and return to the depot; Constraints (13-14) are the TSP constraints for ensuring the connectivity of the tour and preventing sub-tours, respectively; Constraint~(15) is the Knapsack constraint that prevents the vehicle from visiting all the locations.

To apply MSSP to OP, we train linear-SVM to predict whether an edge $(i,j)$ is in the most profitable tour. To initialize the statistical features, random solutions are generated using PS (Algorithm~\ref{alg:ps}) by setting a uniform probability (e.g., $\bm{p}=\bm{1}$) for all decision variables. We employ three problem-specific features. The first one is the distance of the edge normalized by the budget limit, computed by $\frac{d_{i,j}}{T}$. The remaining two features are based on the ratio between the prize gain by visiting a location and the distance of traveling to the location given the current location. The two features are normalized differently, respectively defined as
\begin{equation}
    \frac{p_j / d_{i,j}}{\max_{k=1,\cdots,n} p_k / d_{i,k}},\;\;\;\frac{p_j / d_{i,j}}{\max_{k=1,\cdots,n} p_k / d_{k,j}}.
\end{equation}

Given an unseen OP, the trained linear-SVM can predict whether the edges are likely in the most profitable tour, and the ML prediction is used to guide PS. Similar to TSP, PS starts from the depot location and iteratively samples the next location from a candidate set to visit. Given the current partial tour, the candidate locations are those 1) that are not visited before, and 2) that meet the condition: the total distance from the current position to the candidate vertex and then to the depot is no larger than the remaining budget. PS returns to the depot location to form a feasible OP tour when the candidate set is empty.

\subsection{Maximum Weight Clique Problem}
\label{subsec:mwcp}
MWCP aims at finding a clique on a graph such that the sum of vertex weights of that clique is maximized. Let $G (\mathcal{V}, \mathcal{E}, \mathcal{W})$ be an un-directed weighted graph, consisting of a set of vertices $i \in \mathcal{V}$, a set of edges $(i, j)\in \mathcal{E}$, and a set of weights $w_i \in \mathcal{W}$ each of which corresponding to a vertex $i$. Let $\overline{\mathcal{E}} = \{(i,j)\,|\,\forall (i,j)\not\in \mathcal{E}\}$ denote the set of edges not in $\mathcal{E}$ and let the decision variable $x_i \in \{0,1\}$ indicate whether a vertex $i$ is a part of the clique. The MWCP can be defined as

\begin{align}
    \max_{\bm{x}} \; & \sum_{i \in \mathcal{V}} w_i x_i, \\
    s.t.\; & x_i + x_j \le 1;\; && (i,j) \in \overline{\mathcal{E}}, \\
          & x_i \in \{0,1\};\; && i=\{1,2,\cdots, |\mathcal{V}|\}.
\end{align}

To apply MSSP to MWCP, we train linear-SVM to predict whether a vertex $i$ is a part of the maximum weight clique. To initialize the statistical features, we generate random cliques using PS with uniform probabilities. In addition, we adopt several features to characterize a vertex~$i$, including the weight of the vertex $w_i$; the degree of the vertex, defined by the number of its neighbors $|\mathcal{N}_i|$; the upper bound of the clique containing the vertex $i$, computed by $w_i + \sum_{j \in \mathcal{N}_i} w_j$; the graph density, computed by $2|\mathcal{E}|/(|\mathcal{V}|(|\mathcal{V}|-1))$.

For an unseen MWCP, the trained linear-SVM can predict the vertices likely to be in the maximum weight clique. The ML prediction is used to guide PS to construct a clique by sampling promising vertices from the candidate set~$\mathcal{C}$. The vertices in the candidate set must be neighbors to the vertices already in the partial clique, and the candidate set is constructed by an efficient algorithm based on the set representation of cliques~\citep{sun2019using}.

\subsection{Problem Instances and Training}
\label{subsec:mssp-param}

We generate problem instances for each COP following related studies. For TSP, we generate its instances by sampling a set of locations uniformly randomly in a unit square~\citep{bello2016neural,joshi2019efficient, kool2018attention, kobeaga2018efficient}. For OP, we generate a set of locations the same way as TSP; the distance budget~$T$ is uniformly sampled from $[\frac{d}{2}-1, \frac{d}{2}+1]$ to maximize the problem difficulty~\citep{fischetti1998solving}, where $d$ is the average objective value of optimal solutions for solving TSP of a particular size; the prize attached to a location can be generated by a certain scheme chosen from three schemes, constant, uniform, and distance~\citep{fischetti1998solving}. For MWCP, we generate graphs with a certain size and density using the well-known Barabási-Albert (BA) method~\citep{barabasi1999emergence} and Erdős-Rényi (ER) method~\citep{erdHos1960evolution}, and we assign weights to graph vertices according to~\citep{jiang2018two,cai2016fast}.

\begin{table}[th!]
\centering
\caption{The type of decision variables and the problem-specific features in linear-SVM for three tested problems.}
\label{tab:f}
\resizebox{0.75\columnwidth}{!}{
\begin{tabular}{@{}ccc@{}}
\toprule
Problem & Decision variable & Problem-specific features\\
\midrule
TSP & edge & normalized edge distances\\ \midrule
OP & edge & normalized edge distance, normalized per-cost gains \\ \midrule
MWCP & vertex & graph density, vertex weight, vertex degree, vertex upper bound \\
\bottomrule
\end{tabular}}
\end{table}

For a test problem, we solve $100$ small problem instances and use them to construct the training data. Note that each training example corresponds to a decision variable, consisting of its optimal solution value as the label and its features including statistical features and problem-specific ones (Table~\ref{tab:f}). We train linear-SVM following the standard procedure~\citep{chang2011libsvm}. We raise the value of the penalty parameter for misclassifying positive examples~$r^{+}$ to the inverse of the ratio between the number of positive and negative examples in the training data, whereas the penalty parameter for misclassifying negative examples~$r^{-}$ is set to~$1$. This is because training data can be highly imbalanced, e.g., for TSP the number of edges in the optimal tour is much smaller than the number of non-optimal edges. Adjusting the penalty parameters in such a way encourages linear-SVM to correctly classify positive training examples, which in our situation is more critical than correctly classifying negative ones~\citep{quesnelWDS22}. For the sample size in PS, we set it as a function of the problem size empirically. That is $20n$ for TSP, $50n$ for OP, and $n$ for MWCP, where $n$ denotes the number of vertices (or cities) in a graph. We use Bayesian Optimization (BO) \citep{snoek2012practical} for tuning the parameters~$\beta$ in the logistic model. The BO software~\citep{fernando2014bo} runs MSSP with different sets of~$\beta$ values over~$30$ training instances and identifies the $\beta$ parameters that yield the best solution quality.

\section{Computational Results}
\label{sec:ret}

In this section, we demonstrate that 1) MSSP can effectively improve the ML prediction so as to guide the search to find better-quality solutions (Section~\ref{subsec:mssp-exp}), and 2) MSSP is very competitive compared to ML-based methods (Section~\ref{subsec:exp-ml}) and other optimization methods~(Section~\ref{subsec:exp-trad}).

To assess MSSP and compare it with other methods, we report the gap between the best found solution of a method and the optimal solution of a problem instance, defined as
\begin{equation}
\text{primal gap} \coloneqq \left|\frac{obj^{best\_found}}{obj^{optimal}} - 1\right| \times 100\%.
\end{equation}
In addition, we use average precision~\citep{zhu2004recall} to assess the quality of ML predictions for solution-prediction-based methods following previous work~\citep{ding2020accelerating,shen21learning}. Average precision is computed as the area under a curve that measures the trade-off between precision and recall at different decision thresholds. We note that both measures (i.e., primal gap and the average precision) are computed using the optimal solution of the test problem instance. When the optimal solution is unknown, we report the objective value of the best found solution directly.

\subsection{Multi-shot Solution Prediction}
\label{subsec:mssp-exp}

\begin{figure}[h!]
    \centering
	\begin{tikzpicture}
	\begin{groupplot}[group style = {group size = 3 by 1, horizontal sep = 40pt}, height=0.35\textwidth, width=0.35\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel =  \footnotesize $k^{th}$ iteration of MSSP, axis y line*=left, ylabel =  \footnotesize Average precision of the ML prediction, y label style={at={(axis description cs:-0.12,.5)},anchor=south},x label style={at={(axis description cs:0.5,-0.12)},anchor=north}, xtick = {1,5,10,15,20}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={at={(0.6,0.1)},anchor=south, column sep = 1pt,font=\scriptsize, legend columns = 1,draw=none}, title style={yshift=-1.5ex,}]
   \nextgroupplot[%
    title=\footnotesize {MWCP (ER graph, $1000$ nodes, density $0.18$)},
    ]
	\addplot[color=red, mark=square*, mark size = 1, line width=0.30mm] table [x=x, y=ER_0.18, col sep=comma] {figures/mwc_ap.txt};
   \nextgroupplot[%
    title=\footnotesize {TSP ($100$ cities)},
    ]
	\addplot[color=red, mark=square*, mark size = 1, line width=0.30mm] table [x=x, y=y100, col sep=comma] {figures/tsp_ap.txt};
	\nextgroupplot[%
    title=\footnotesize OP ($50$ locations and `uniform' price),
    ]
	\addplot[color=red, mark=square*, mark size = 1, line width=0.30mm] table [x=x, y=y50unif, col sep=comma] {figures/op_ap.txt};
    \end{groupplot}
	\end{tikzpicture}
\caption{The average precision of the ML prediction at every iteration of MSSP for three tested COPs. For each COP, the mean statistics across $100$ test instances are reported, and information for the test instances is shown in brackets. In MSSP, the linear-SVM and the logistic model are trained using the solved problem instances of the same size as the test ones. }
\label{fig:mssp-ap}
\end{figure}

\begin{figure}[h!]
\vspace{3.em}
	\begin{minipage}[c][1\width]{
	   0.45\textwidth}
	   \centering
	\begin{tikzpicture}
	\begin{groupplot}[group style = {group size = 1 by 1, horizontal sep = 50pt}, height=\textwidth,width=\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel = \scriptsize $k^{th}$ iteration of MSSP , xtick = {1,5,10,15,20},x label style={at={(axis description cs:0.5,-0.06)},anchor=north}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={column sep = 1pt, legend columns = 1,font=\scriptsize,draw=none}]
   \nextgroupplot[%
    ,legend to name=grouplegend,
    title= MWCP,
    ylabel = \scriptsize Average precision of the ML prediction,
    y label style={at={(axis description cs:-0.15,.5)},anchor=south},
    ]
	\addplot[color=orange, mark=square, mark size = 1, line width=0.30mm] table [x=x, y=ER_0.1, col sep=comma] {figures/mwc_ap.txt}; \addlegendentry{Density: $0.10$}
	\addplot[color=blue, mark=triangle, mark size = 1, line width=0.30mm] table [x=x, y=ER_0.18, col sep=comma] {figures/mwc_ap.txt}; \addlegendentry{Density: $0.18$}
	\addplot[color=red, mark=diamond, mark size = 1, line width=0.30mm] table [x=x, y=ER_0.26, col sep=comma] {figures/mwc_ap.txt}; \addlegendentry{Density: $0.26$}
    \end{groupplot}
    \node at (group c1r1.west) [anchor=west, yshift=-1.75cm, xshift=2.75cm] {\ref{grouplegend}};
	\end{tikzpicture}
\caption{On MWCP, the average precision of the ML prediction at every iteration of MSSP. The test problem instances are of three different graph densities, and the mean statistics across $100$ test instances for each one are reported. In MSSP, the linear-SVM and the logistic model are trained using the solved problem instances of $1000$ nodes with a density value~$0.18$.}

\label{fig:mwcp-generalize}	\end{minipage}
 \hfill
	\begin{minipage}[c][1\width]{
	   0.45\textwidth}
	   \centering
	\begin{tikzpicture}
	\begin{groupplot}[group style = {group size = 1 by 1, horizontal sep = 50pt}, height=\textwidth,width=\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel = \scriptsize $k^{th}$ iteration of MSSP, xtick = {1,5,10,15,20},x label style={at={(axis description cs:0.5,-0.06)},anchor=north}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={column sep = 1pt, legend columns = 1,font=\scriptsize,draw=none}]
   \nextgroupplot[%
    ,legend to name=grouplegend,
    title=TSP,
    ylabel = \scriptsize Average precision of the ML prediction,
    y label style={at={(axis description cs:-0.15,.5)},anchor=south},
    ]
	\addplot[color=blue, mark=triangle, mark size = 1, line width=0.30mm] table [x=x, y=y100, col sep=comma] {figures/tsp_ap.txt}; \addlegendentry{\# cities: $100$}
	\addplot[color=marine, mark=square*, mark size = 1, line width=0.30mm] table [x=x, y=y300, col sep=comma] {figures/tsp_ap.txt}; \addlegendentry{\# cities: $300$}
	\addplot[color=red, mark=*, mark size = 1, line width=0.30mm] table [x=x, y=y1000, col sep=comma] {figures/tsp_ap.txt}; \addlegendentry{\# cities: $1000$}
    \end{groupplot}
    \node at (group c1r1.west) [anchor=west, yshift=-1.75cm, xshift=2.75cm] {\ref{grouplegend}};

	\end{tikzpicture}
\caption{ On TSP, the average precision of the ML prediction at every iteration of MSSP. The test problem instances are of three different sizes, and the mean statistics across $100$ test instances for each size are reported. In MSSP, the linear-SVM and the logistic model are trained using the solved problem instances of $100$ cities.}

\label{fig:tsp-generalize}	\end{minipage}
\vspace{3.em}
\end{figure}

\def\paragraph#1{\noindent\emph{#1} }
\paragraph{Can MSSP improve the quality of ML predictions?} Figure~\ref{fig:mssp-ap} shows the average precision of the ML predictions over iterations of MSSP. For each COP, MSSP is trained using a certain size of problem instances and is tested on $100$ instances of the same size. As can be seen, on all the three tested problems MSSP can effectively improve the quality of the ML predictions particularly at the first several iterations. Note that the improvement in prediction quality can be attributed to MSSP's superb ability to iteratively refine statistical features using high-quality solutions, i.e., MSSP improves the quality of prediction via harnessing results from the search. In addition, MSSP can also effectively improve the ML prediction of linear-SVM when the training instances and test instances are of different sizes (but still the same problem), as shown in Figure \ref{fig:mwcp-generalize} for MWCP and Figure \ref{fig:tsp-generalize} for TSP.

\begin{figure}[h!]
    \centering
	\begin{tikzpicture}
	\begin{groupplot}[group style = {group size = 3 by 1, horizontal sep = 40pt}, height=0.35\textwidth, width=0.35\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel =  \footnotesize $k^{th}$ iteration of MSSP, axis y line*=left, ylabel =  \footnotesize Primal gap, y label style={at={(axis description cs:-0.12,.5)},anchor=south},x label style={at={(axis description cs:0.5,-0.12)},anchor=north}, xtick = {1,5,10,15,20}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={at={(0.8,0.8)},anchor=south, column sep = 1pt,font=\scriptsize, legend columns = 1,draw=none}, title style={yshift=-1.5ex,}]
   \nextgroupplot[%
    title=\footnotesize MWCP (ER graph with density $0.18$),
    ]
	\addplot[color=marine, mark=triangle*, mark size = 1.5, line width=0.30mm] table [x=x, y=ER_0.18, col sep=comma] {figures/mwc_gap.txt};
   \nextgroupplot[%
    title=\footnotesize TSP ($100$ cities),
    ]
	\addplot[color=marine, mark=triangle*, mark size = 1.5, line width=0.30mm] table [x=x, y=y100, col sep=comma] {figures/tsp_greedy_obj.txt};

	\nextgroupplot[%
    title=\footnotesize OP ($50$ locations and `uniform' price),
    ]
    \addplot[color=marine, mark=triangle*, mark size = 1.5, line width=0.30mm] table [x=x, y=y50unif, col sep=comma] {figures/op_greedy_obj_opt.txt};
    \end{groupplot}
	\end{tikzpicture}
\caption{The primal gap of the best solution at every iteration of MSSP for three tested COPs. For each COP, the mean statistics across $100$ test instances are reported, and the size of the test problem instances is shown in brackets. In MSSP, the linear-SVM and the logistic model are trained using the solved problem instances of the same size as the test ones. Note that the primal-gap curve is not guaranteed to monotonically decrease, as we report the best solution constructed using the ML prediction at a particular iteration rather than the best found solution so far.}
\label{fig:mssp-gap}
\end{figure}

\paragraph{Can improved ML predictions guide the search to find better solutions? }
Figure~\ref{fig:mssp-gap} shows the primal gap of the best solution at each iteration of MSSP. We can observe that the search method can construct much better solutions as a result of making use of the iteratively improved ML prediction. Furthermore, the improvements in the solution quality are particularly large at the first several iterations. This observation is consistent with the fast improvement of the prediction quality as shown in Figure~\ref{fig:mssp-ap}. These results show that the effectiveness of the ML prediction and the quality of the best found solution are closely related, and improving the quality of ML prediction is critical in boosting a search method.

\begin{figure}[h!]
    \centering
	\begin{tikzpicture}
	\begin{axis} [legend style={nodes={scale=0.5, transform shape}}, title=  MWCP, box plot width=0.0em,  height=0.33\textwidth,width=0.33\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel = \scriptsize Graph density, ylabel = \scriptsize Milliseconds, y label style={at={(axis description cs:-0.12,.5)},anchor=south}, xtick = {0.1, 0.18, 0.26},x label style={at={(axis description cs:0.5,-0.12)},anchor=north}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={at={(0.47,0.22)},anchor=west,font=,draw=none},
	mark repeat=1,mark phase=0]
	\addplot[color=black, line width=0.30mm, error bars/.cd, y dir=both, y explicit] table [x=name, y=mean, y error=std, col sep=comma] {figures/mwc_time.txt};
	\end{axis}
	\end{tikzpicture}
	\hspace{2em}
	\begin{tikzpicture}
	\begin{axis} [box plot width=0.0em, title=  TSP,  height=0.33\textwidth,width=0.33\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel = \scriptsize Graph size, ylabel = \scriptsize Seconds, y label style={at={(axis description cs:-0.12,.5)},anchor=south}, xtick = {100,500, 1000},x label style={at={(axis description cs:0.5,-0.12)},anchor=north}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={at={(0.6,0.24)},anchor=west,font=,draw=none},
	mark repeat=1,mark phase=0]
	\addplot[color=black, line width=0.30mm, error bars/.cd, y dir=both, y explicit] table [x=name, y=mean, y error=std, col sep=comma] {figures/tsp_time.txt};
	\end{axis}
	\end{tikzpicture}
	\hspace{2em}
	\begin{tikzpicture}
	\begin{axis} [box plot width=0.0em,  height=0.33\textwidth,width=0.33\textwidth, title=  OP, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel = \scriptsize Graph size, ylabel = \scriptsize Seconds, y label style={at={(axis description cs:-0.12,.5)},anchor=south}, xtick = {100,300,500},x label style={at={(axis description cs:0.5,-0.12)},anchor=north}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={at={(0.6,0.24)},anchor=west,font=,draw=none},
	mark repeat=1,mark phase=0]
	\addplot[color=black, line width=0.30mm, error bars/.cd, y dir=both, y explicit] table [x=name, y=mean, y error=std, col sep=comma] {figures/op_time.txt};
	\end{axis}
	\end{tikzpicture}
\caption{Time spent per iteration when problem size grows.}
\label{fig:time}
\end{figure}

\paragraph{Computational time for MSSP.} Figure~\ref{fig:time} shows the computational time spent per MSSP iteration when problem size grows. On MWCP, MSSP only spends a few milliseconds per iteration and its computational time increases slowly with respect to the number of edges in a graph. For TSP and OP, the computational time per MSSP iteration is in seconds and increases in polynomial time with respect to the number of cities. These observations can be explained by analyzing the PS method, which takes the majority of the computational time of MSSP in practice. The time complexity of PS is determined by the sample size, the solution length, and the time complexity for resolving conflict variables. For MWCP, its sample size and the length of a solution are much smaller than those for TSP and OP. When problem size increases, the solution length for MWCP grows much slower than that for TSP and OP.

\subsection{Comparison to learning-based methods}
\label{subsec:exp-ml}

This part compares MSSP with several learning-based methods, including: 1) SSSP-SVM, a special case of MSSP with the number of iterations~$T=1$. More specifically, the components of SSSP-SVM are identical to MSSP's, but SSSP-SVM only predicts the optimal solution once and uses the ML prediction to guide the search. 2) SSSP-GNN, which refers to the method developed for TSP\footnote{Code: \url{https://github.com/chaitjo/graph-convnet-tsp}.}~\citep{joshi2019learning} and the one for the maximum clique problem\footnote{Code: \url{https://github.com/intel-isl/NPHard}.}~\citep{li2018combinatorial}. Note that these methods both use the popular Graph Neural Network~(GNN) to predict optimal solutions, but they are different in the use of (problem-specific) search methods. We also note that \citet{li2018combinatorial} initially tested their approach on a special case of MWCP with the uniform vertex weight~$w_i=1$; hence their feature vectors for vertices are one-hot $\bm{x}=\bm{1}$. In our case of MWCP, vertices have different weights; hence the feature vector of a vertex~$i$ is set to~$\bm{x_i}=\bm{w_i}$ accordingly. 3) RL-Attn\footnote{Code: \url{https://github.com/wouterkool/attention-learn-to-route}.}, a RL-based approach designed for routing problems~\citep{kool2018attention}. RL-Attn employs an attention-based neural network (Attn) as the policy, which is trained end-to-end via RL to construct high-quality solutions directly.

For testing, a method trained on a COP is used to solve unseen problem instances of the same COP. We parallelize the sampling process of our MSSP and SSSP-SVM on $8$ CPUs (Intel(R) Core(TM) i5-8300H CPU @ 2.30GHz). For deep-learning methods, we equip them with a Tesla P100 GPU, and their computational time can be substantially reduced as compared to that when using $8$ CPUs.

\begin{table*}[h!]
\centering
\caption{The results of the learning-based methods on MWCP. For a type of graph, the training set is formed from solved instances with $1000$ vertices and a density value of $0.18$, and the test instances are formed on the same type of graph but of different density values. A method is set to generate a total number of $10n$ solution. The primal gap of the best solution is reported, averaged over $100$ test instances. The best results are highlighted in bold (statistically significant based on a $t$-test, with the significance level set to $0.05$).}
\label{tab:mwc}
\resizebox{0.95\textwidth}{!}{
\begin{tabular}{@{}cccccccccc@{}}
\toprule
\multirow{2}{*}{Graph type} & \multirow{2}{*}{Density} & \multicolumn{3}{c}{Average precision} & \multicolumn{3}{c}{Primal gap} & \multicolumn{2}{c}{Time} \\
& & MSSP & SSSP-SVM & SSSP-GNN & MSSP & SSSP-SVM & SSSP-GNN & MSSP/SSSP-SVM & SSSP-GNN \\
\cmidrule(lr){1-2}\cmidrule(lr){3-5}\cmidrule(lr){6-8}\cmidrule(lr){9-10}
\multirow{3}{*}{Erdős-Rényi}& $0.10$ & \textbf{0.17} & 0.08 & 0.09 & 0.0 & 0.04 & 1.19 & 0.1 & 0.1\\
& $0.18$ & \textbf{0.23} & 0.14 & 0.10 & 0.0 & 0.0 & 3.74 & 0.1 & 0.1 \\
& $0.26$  & \textbf{0.21} & 0.12 & 0.09 & \textbf{0.17} & 0.23 & 6.71 & 0.1 & 0.2\\
\cmidrule(lr){1-2}\cmidrule(lr){3-5}\cmidrule(lr){6-8}\cmidrule(lr){9-10}
\multirow{3}{*}{Barabási-Albert}& $0.10$ & \textbf{0.95} & 0.84 & 0.49 & 0.0 & 0.0 & 3.78 & 0.1 & 0.2\\
& $0.18$  & \textbf{0.98} & 0.93 & 0.88 & 0.0 & 0.03 & 0.14 & 0.4 & 0.5\\
& $0.26$ & \textbf{0.99} & 0.97 & 0.53 & 0.0 & 0.06 & 8.82 & 0.2 & 0.3 \\
\bottomrule
\end{tabular}}
\end{table*}

Table~\ref{tab:mwc} shows the results on MWCP. It can be seen that MSSP produces the best-quality ML prediction among the three solution-prediction-based methods. Further, we observe the performances of both MSSP and SSSP are robust, i.e., their performance does not significantly degrade when testing on the graph with a different density to the training ones. Being consistent with the average precision, the best found solution by MSSP has better quality than that of other methods.

\begin{table*}[h!]
\centering
\caption{The results of the learning-based methods on TSP. All methods are trained using TSP instances of $100$ cities, and are tested on TSP instances of different sizes. The cutoff time for MSSP and SSSP-SVM is set to the minimum computational time between SSSP-GNN and RL-Attn under their default parameter settings. The primal gap of the best solution is reported, averaged over $100$ test instances. The best results are highlighted in bold. }
\label{tab:tsp}
\resizebox{0.95\textwidth}{!}{
\begin{tabular}{@{}cccccccccccc@{}}
\toprule
\multirow{2}{*}{Graph size}& \multicolumn{3}{c}{Average precision}& \multicolumn{4}{c}{Primal gap} & \multicolumn{3}{c}{Time} \\
& MSSP & SSSP-SVM & SSSP-GNN & MSSP & SSSP-SVM & SSSP-GNN & RL-Attn & MSSP/SSSP-SVM & SSSP-GNN & RL-Attn \\
\cmidrule(lr){1-1} \cmidrule(lr){2-4}\cmidrule(lr){5-8}\cmidrule(lr){9-11}
$30$ & \textbf{0.84} & 0.75& 0.71 & \textbf{0.6} & 1.0 & 33.6 & 2.2 & 0.1 & 0.1 & 0.1 \\
$100$ & 0.85 & 0.65& \textbf{0.98} & 4.1 & 11.9 & \textbf{2.1}  & \textbf{2.3} & 0.4 & 0.4 & 0.7 \\
$300$ & \textbf{0.85} & 0.52 & 0.42 & 12.0 & 32.9 & 71.5 & \textbf{10.4} & 3.8 & 3.8 &5.8 \\
$500$ & \textbf{0.84} & 0.46 & 0.34 & \textbf{14.1} & 41.5& 79.1 & 17.6 & 12.4 & 12.4 & 16.7 \\
$1000$ & \textbf{0.83} & 0.36 & 0.23 & \textbf{16.8} & 55.8 & 117.7 & 29.0 &48.6 &48.6 & 75.7 \\
\bottomrule
\end{tabular}}
\end{table*}

Table~\ref{tab:tsp} presents the results on TSP, showing that MSSP scales well across TSP instances of different sizes as evidenced by the average precision, while SSSP-SVM performs reasonably well only for small-scale TSP instances with 30 cities. Consequently, the quality of the best solution found by MSSP is much more close to the optimal one than that found by SSSP-SVM. It is interesting to note that for SSSP-GNN and RL-Attn, their performances are competitive only when the size of the training instances (i.e., $100$ cities) and testing instances are similar. However, their performances deteriorate quickly when testing on larger TSP instances.

\begin{table*}[h!]
\centering
\caption{The results of the learning-based methods on OP. All methods are trained using OP instances of $50$ locations and the uniform prize, and are tested on OP instances of different sizes and different prize types. The cutoff time for MSSP and SSSP-SVM is set to the average computational time of RL-Attn under its default parameter setting. OP is a maximization problem, and the objective value of the best solution is reported, averaged over $100$ test instances. The best results are highlighted in bold.}
\label{tab:op}
\resizebox{0.95\textwidth}{!}{
\begin{tabular}{@{}cccccccccccc@{}}
\toprule
\multirow{2}{*}{Graph size} & \multicolumn{3}{c}{Prize type: `uniform'} & \multicolumn{3}{c}{Prize type: `constant'} & \multicolumn{3}{c}{Prize type: `distance'} & \multicolumn{2}{c}{Time}\\
& MSSP& SSSP-SVM& RL-Attn & MSSP & SSSP-SVM & RL-Attn & MSSP & SSSP-SVM &RL-Attn& MSSP/SSSP-SVM & RL-Attn\\
\cmidrule(lr){1-1} \cmidrule(lr){2-4}\cmidrule(lr){5-7}\cmidrule(lr){8-10}\cmidrule(lr){11-12}

$50$  & 14.9 & 14.8 & 15.0 & 27.1 & 26.9 & 26.5 & 14.9 & 14.7 & 14.9 & 0.2 & 0.2 \\
$100$ & 30.8 & 29.2 & 31.2 & \textbf{56.4} & 53.8 & 54.9 & 30.8 & 29.2 & 30.7 & 0.9 & 0.9\\
$300$ & 79.4 & 70.7 & \textbf{81.6} & \textbf{147.3} & 133.9 & 143.3 & \textbf{83.1} & 72.5 & 79.0 & 6.5 & 6.5\\
$500$ & \textbf{125.0} & 112.8 & 111.9 & \textbf{232.7} & 218.0 & 201.2 & \textbf{131.5} & 105.3 & 100.9 & 14.4 & 14.4\\

\bottomrule
\end{tabular}}
\end{table*}

Table~\ref{tab:op} shows the results on OP. Overall, MSSP significantly outperforms other methods in most cases, though RL-Attn wins with statistical significance on just one occasion. Similar to the observations on TSP, MSSP significantly improves on SSSP-SVM especially for large OP instances. Considering that the underlining distribution of the real-world problems is usually unknown, these results further demonstrate that MSSP is more generally applicable to real-world scenarios.

\subsection{Comparison to traditional methods}
\label{subsec:exp-trad}

We mainly compare MSSP with ant colony system~(ACS)~\citep{dorigo1996ant}, a representative in the class of EDAs~(Section~\ref{subsec:MHcompare}). ACS has been extensively studied in combinatorial optimization, and we adopt ACS variants in ~\citep{xu2007improved, ke2008ants, dorigo1997ant} for our test problems. In addition to ACS, we include Gurobi \citep{gurobi2018gurobi} as a baseline, a commercial mixed-integer-programming solver with a bag of primal heuristics. As our focus is on solution quality, we instruct Gurobi to spend $90\%$ of the computational time on executing these primal heuristics by setting the `Heuristics' parameter to $0.9$.

\begin{figure*}[t!]
    \centering
	\begin{tikzpicture}
	\begin{groupplot}[group style = {group size = 3 by 2, horizontal sep = 45pt, vertical sep = 45pt}, height=0.3\textwidth, width=0.3\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel = \scriptsize Seconds, ylabel = \scriptsize Primal gap, y label style={at={(axis description cs:-0.1,.5)},anchor=south},x label style={at={(axis description cs:0.5,-0.05)},anchor=north}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={column sep = 1pt,font=\scriptsize, legend columns = 1,draw=none}, title style={yshift=-1.5ex,}]
   \nextgroupplot[%
    legend to name=grouplegend1,
    title=\footnotesize 30 cities,
    ]
	\addplot[color=red, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/mssp_30.txt}; \addlegendentry{MSSP}
	\addplot[color=blue, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/aco_30.txt}; \addlegendentry{ACS}
	\addplot[color=black, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/gurobi_30.txt}; \addlegendentry{Gurobi}
   \nextgroupplot[%
    title=\footnotesize100 cities,
    ]
	\addplot[color=red, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/mssp_100.txt};
	\addplot[color=blue, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/aco_100.txt};
	\addplot[color=black, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/gurobi_100.txt};
	\nextgroupplot[%
    title=\footnotesize300 cities,
    ]
	\addplot[color=red, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/mssp_300.txt};
	\addplot[color=blue, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/aco_300.txt};
	\addplot[color=black, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/gurobi_300.txt};
	\nextgroupplot[%
    title=\footnotesize500 cities,
    ]
	\addplot[color=red, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/mssp_500.txt};
	\addplot[color=blue, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/aco_500.txt};
	\addplot[color=black, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/gurobi_500.txt};
	\nextgroupplot[%
    title=\footnotesize1000 cities,
    ]
	\addplot[color=red, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/mssp_1000.txt};
	\addplot[color=blue, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/aco_1000.txt};
	\addplot[color=black, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/tsp/gurobi_1000.txt};
    \end{groupplot}
    \node at (group c3r1.south) [anchor=south, yshift=-3.5cm, xshift=0cm] {\ref{grouplegend1}};
	\end{tikzpicture}
\caption{The results on TSP using MSSP, ACS and Gurobi. For each size of TSP, the results are averaged over $100$ test instances.}
\label{fig:tsp-extended}
\end{figure*}
\begin{figure*}[t!]
    \centering
	\begin{tikzpicture}
	\begin{groupplot}[group style = {group size = 4 by 1, horizontal sep = 45pt, vertical sep = 45pt}, height=0.25\textwidth, width=0.25\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel = \scriptsize Seconds, ylabel = \scriptsize Objective, y label style={at={(axis description cs:-0.12,.5)},anchor=south},x label style={at={(axis description cs:0.5,-0.05)},anchor=north}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={column sep = 1pt,font=\scriptsize, legend columns = 4,draw=none},  title style={yshift=-1.5ex}]
   \nextgroupplot[%
    ,legend to name=grouplegend,
    title=\footnotesize50 locations,
    ]
	\addplot[color=red, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/mssp_50unif.txt}; \addlegendentry{MSSP}
	\addplot[color=blue, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/aco_50unif.txt}; \addlegendentry{ACS}
	\addplot[color=black, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/gurobi_50unif.txt}; \addlegendentry{Gurobi}
   \nextgroupplot[%
    title=\footnotesize100 locations,
    ]
	\addplot[color=red, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/mssp_100unif.txt};
	\addplot[color=blue, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/aco_100unif.txt};
	\addplot[color=black, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/gurobi_100unif.txt};
	\nextgroupplot[%
    title=\footnotesize300 locations,
    ]
	\addplot[color=red, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/mssp_300unif.txt};
	\addplot[color=blue, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/aco_300unif.txt};
	\addplot[color=black, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/gurobi_300unif.txt};
	\nextgroupplot[%
    title=\footnotesize500 locations,
    ]
	\addplot[color=red, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/mssp_500unif.txt};
	\addplot[color=blue, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/aco_500unif.txt};
	\addplot[color=black, line width=0.30mm] table [x=x, y=y, col sep=comma] {figures/op/gurobi_500unif.txt};

    \end{groupplot}

    \node at (group c1r1.north) [anchor=west, yshift=0.75cm, xshift=3.5cm] {\ref{grouplegend}};
	\end{tikzpicture}
\caption{The results on OP~(Maximization) using MSSP, ACS and Gurobi. The test OP instances are of prize type `uniform'. For each size of OP, the results are averaged over $100$ test instances. For MSSP, the logistic model in MSSP is optimized using OP instances with $500$ locations. }
\label{fig:op-extended}
\end{figure*}

The results for the two routing problems, TSP and OP, are shown in Figure~\ref{fig:tsp-extended} and \ref{fig:op-extended}, respectively. For both test problems, MSSP can find high-quality solutions more quickly, within just a few seconds. In contrast, ACS tends to converge to solutions that are worse than that of MSSP, and Gurobi is very effective in tackling small TSP instances but does not work well for larger ones.

\begin{figure*}[t!]
    \centering
	\begin{tikzpicture}
	\begin{groupplot}[group style = {group size = 3 by 3, horizontal sep = 50pt, vertical sep = 50pt}, height=0.32\textwidth, width=0.32\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel = \scriptsize Seconds, ylabel = \scriptsize Objective, y label style={at={(axis description cs:-0.16,.5)},anchor=south}, xtick = {1,5,10,15,20},x label style={at={(axis description cs:0.5,-0.05)},anchor=north}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={column sep = 1pt,font=\scriptsize, legend columns = 1,draw=none},title style={yshift=-1.5ex}]
   \nextgroupplot[%
   legend to name=grouplegend1,
    title= \footnotesize C1000-9,
    ]
	\addplot[color=red, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C1000-9/mssp.txt}; \addlegendentry{MSSP}
	\addplot[color=purple, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C1000-9/sssp.txt}; \addlegendentry{SSSP-SVM}
	\addplot[color=blue, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C1000-9/aco.txt}; \addlegendentry{ACS}
	\addplot[color=gray, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C1000-9/gurobi.txt}; \addlegendentry{Gurobi}
	\addplot[color=black, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/C1000-9/lscc.txt};  \addlegendentry{LSCC+BMS}
	\addplot[color=orange, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/C1000-9/fastwclq.txt};  \addlegendentry{FastWClq}
	\addplot[color=green, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/C1000-9/tsm.txt};\addlegendentry{TSM}
   \nextgroupplot[%
    title= \footnotesize C2000-5,
    ]
	\addplot[color=red, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C2000-5/mssp.txt};
	\addplot[color=purple, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C2000-5/sssp.txt};
	\addplot[color=blue, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C2000-5/aco.txt};
	\addplot[color=gray, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C2000-5/gurobi.txt};
	\addplot[color=black, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma, dashed] {dimacs/C2000-5/lscc.txt};
	\addplot[color=orange, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/C2000-5/fastwclq.txt};
	\addplot[color=green, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/C2000-5/tsm.txt};
   \nextgroupplot[%
    title= \footnotesize C2000-9,
    ]
	\addplot[color=red, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C2000-9/mssp.txt};
	\addplot[color=purple, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C2000-9/sssp.txt};
	\addplot[color=blue, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C2000-9/aco.txt};
	\addplot[color=gray, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C2000-9/gurobi.txt};
	\addplot[color=black, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/C2000-9/lscc.txt};
	\addplot[color=orange, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/C2000-9/fastwclq.txt};
	\addplot[color=green, line width=0.50mm,dashed] table [x=x, y=y, col sep=comma] {dimacs/C2000-9/tsm.txt};
   \nextgroupplot[%
    title= \footnotesize C4000-5,
    ]
	\addplot[color=red, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C4000-5/mssp.txt};
	\addplot[color=purple, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C4000-5/sssp.txt};
	\addplot[color=blue, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C4000-5/aco.txt};
	\addplot[color=gray, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/C4000-5/gurobi.txt};
	\addplot[color=black, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/C4000-5/lscc.txt};
	\addplot[color=orange, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/C4000-5/fastwclq.txt};
	\addplot[color=green, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/C4000-5/tsm.txt};
   \nextgroupplot[%
    title= \footnotesize hamming10-4,
    ]
	\addplot[color=red, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/hamming10-4/mssp.txt};
	\addplot[color=purple, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/hamming10-4/sssp.txt};
	\addplot[color=blue, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/hamming10-4/aco.txt};
	\addplot[color=gray, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/hamming10-4/gurobi.txt};
	\addplot[color=black, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/hamming10-4/lscc.txt};
	\addplot[color=orange, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/hamming10-4/fastwclq.txt};
	\addplot[color=green, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/hamming10-4/tsm.txt};
   \nextgroupplot[%
    title= \footnotesize keller6,
    ]
	\addplot[color=red, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/keller6/mssp.txt};
	\addplot[color=purple, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/keller6/sssp.txt};
	\addplot[color=blue, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/keller6/aco.txt};
	\addplot[color=gray, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/keller6/gurobi.txt};
	\addplot[color=black, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/keller6/lscc.txt};
	\addplot[color=orange, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/keller6/fastwclq.txt};
	\addplot[color=green, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/keller6/tsm.txt};
   \nextgroupplot[%
    title= \footnotesize MANN\_a45,
    ]
	\addplot[color=red, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/MANN-a45/mssp.txt};
	\addplot[color=purple, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/MANN-a45/sssp.txt};
	\addplot[color=blue, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/MANN-a45/aco.txt};
	\addplot[color=gray, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/MANN-a45/gurobi.txt};
	\addplot[color=black, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/MANN-a45/lscc.txt};
	\addplot[color=orange, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/MANN-a45/fastwclq.txt};
	\addplot[color=green, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/MANN-a45/tsm.txt};
   \nextgroupplot[%
    title= \footnotesize p\_hat1500-3,
    ]
	\addplot[color=red, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/p-hat1500-3/mssp.txt};
	\addplot[color=purple, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/p-hat1500-3/sssp.txt};
	\addplot[color=blue, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/p-hat1500-3/aco.txt};
	\addplot[color=gray, line width=0.50mm] table [x=x, y=y, col sep=comma] {dimacs/p-hat1500-3/gurobi.txt};
	\addplot[color=black, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/p-hat1500-3/lscc.txt};
	\addplot[color=orange, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/p-hat1500-3/fastwclq.txt};
	\addplot[color=green, line width=0.50mm, dashed] table [x=x, y=y, col sep=comma] {dimacs/p-hat1500-3/tsm.txt};
    \end{groupplot}
    \node at (group c3r2.south) [anchor=south, yshift=-4.5cm, xshift=0cm] {\ref{grouplegend1}};
	\end{tikzpicture}
\caption{The results generated by the generic and specialized methods on MWCP using large graph instances in the DIMACS benchmark. The results of generic methods are in solid lines and those of specialized methods are in dashed lines. For a method, the mean statistic of best found solutions is reported averaged over $10$ independent runs. We note that the MSSP and SSSP-SVM is trained using a set of $26$ easy graphs with known optimal solutions.}
\label{fig:mwc}
\end{figure*}

Next, we present the results for MWCP using $8$ large graphs from the DIMACS benchmark, as shown in Figure~\ref{fig:mwc}. Here, we focus on the comparison between several generic methods (i.e., those in solid lines). It can be seen that MSSP significantly outperforms Gurobi and ACS, with an exception on MANN\_a45 where Gurobi is very competitive. Notably, MSSP improves upon SSSP-SVM by a large margin on those challenging graphs. Moreover, these results suggest that MSSP is very robust, given that the test graphs in the DIMACS benchmark have diverse characteristics.

We also compare MSSP with specialized methods. As our primary goal is to show that MSSP can improve over the existing SSSP-SVM on \emph{general} COPs, here we only use MWCP as our test problem and compare generic methods against three state-of-the-art solution methods (solvers) for MWCP -- TSM~\citep{jiang2018two}, LSCC+BMS~\citep{wang2016two} and Fastwclq~\citep{cai2016fast}. The results for these specialized methods are shown in dashed lines in Figure~\ref{fig:mwc}. We can observe that MSSP achieves better results than TSM and Fastwclq, and slightly worse results than the LSCC+BMS algorithm. These results are interesting because generic methods (such as MSSP) are usually not expected to be competitive to the specialized methods.

\section{MSSP for Column Generation and Branch-and-Price}
\label{sec:cgbp}

This section applies MSSP as a pricing heuristic for CG~\citep{lubbecke2005selected}, to boost the exact B\&P method~\citep{barnhart1998branch} for solving the graph coloring problem (GCP). CG aims to capture the columns in the optimal solution of a linear program~(LP), and these columns are gradually generated by solving a sequence of NP-hard pricing problems. Using MSSP as the pricing heuristic for solving the pricing problems can accelerate CG in capturing optimal LP columns, because MSSP can generate many high-quality columns efficiently. This section is built on our previous work~\citep{shen2022enhancing}, which proposes a SSSP approach as the pricing heuristic for CG and the exact B\&P method for solving GCP.

\subsection{Background}

In GCP, the goal is 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~\citep{malaguti2010survey}. To solve GCP, B\&P recursively decomposes the original GCP (root node) into subproblems (child nodes) and prunes a node safely without further expansion if its lower bound was no better than the current best found solution. For GCP, a tighter lower bound can be obtained by solving the LP relaxation of its set-covering formulation (GCP-SC), defined as using a minimum number of maximal independent sets~(MISs) to cover all the vertices in a graph~\citep{mehrotra1996column}. Let the binary variable~$z_\mathcal{I}$ indicate whether a MIS~$\mathcal{I} \in \mathbb{I}$ is used to cover the vertices in a graph $G(\mathcal{V},\mathcal{E})$, where $\mathbb{I}$ is the set of all the possible MISs; $\mathcal{V}$ denotes the set of vertices; $\mathcal{E}$ denotes the set of edges. GCP-SC can be formulated as
\begin{align}
    \min_{\bm{z}} \;& \sum_{m \in \mathcal{I}} z_\mathcal{I}, & \text{(GCP-SC)}\\
        s.t. \;& \sum_{\mathcal{I} \in \mathbb{I}, i\in \mathcal{I}} z_\mathcal{I} \geq 1; &  i \in \mathcal{V}, \label{eq:at-least-one} \\
               & z_\mathcal{I} \in \{0,1\}; &\mathcal{I} \in \mathbb{I},
\end{align}
\noindent where Constraints (24) specify that, for each vertex $i$, at least one MIS containing that vertex $i \in \mathcal{I}$ must be in use.

Given GCP-SC, its LP relaxation can be obtained by relaxing integer constraints, i.e., replacing Constraints (25) to $z_\mathcal{I} \in [0,1]; \mathcal{I} \in \mathbb{I}$. Solving the LP can yield a tight lower bound that can help B\&P aggressively prune nodes and potentially reduce its computational time. However, solving such a LP itself can be challenging to existing methods (e.g., the simplex methods or the interior point methods~\citep{dantzig2016linear}), because the LP may need an exponential number of variables (or columns) to represent all MISs in a graph.

CG is an iterative method for solving large-scale LPs, based on the fact that an optimal LP solution typically has only a tiny proportion of optimal LP columns, i.e., columns with non-zero values. CG starts by solving a restricted master problem (RMP) that considers a subset of columns (or MISs) in the original LP and gradually includes new MISs that can improve the solution to the current RMP, i.e., those with negative reduced costs (NRC). Searching for NRC MISs can be formulated as an optimization problem (i.e., the pricing problem), with minimizing the reduced cost as the objective function, defined for GCP-SC as
\begin{align}
    \min_{\bm{x}} \;& 1 - \sum_{i \in \mathcal{V}} \pi_i\cdot x_i, & \text{(MWISP)}\\
        s.t. \;& x_i + x_j \leq 1; & (i,j)\in \mathcal{E}, \\
              & x_i \in \{0,1\}; & i \in \mathcal{V},
\end{align}
\noindent where $x_i$ denotes whether a vertex $i$ is a part of the MIS and $\pi_i$ holds the optimal dual solution associated to that vertex $i$ to the current RMP. As can be seen, the pricing problem is the maximum weight independent problem (MWISP) on the graph $G(\mathcal{V},\mathcal{E}, \mathcal{W})$ with optimal dual solution as vertex weights $\mathcal{W} \coloneqq \{\pi_i\;|\;i \in \mathcal{V}\}$, which is equivalent to MWCP on its complementary graph $G(\mathcal{V},\overline{\mathcal{E}}, \mathcal{W})$ as shown in Equations~(19)-(21).

CG alternates between solving the RMP and the MWISP, and the computational process terminates when no NRC column exists. This certifies that the RMP has already captured all the optimal LP columns to the original LP; hence the optimal solution to the RMP is also optimal to the original LP.

\subsection{Experiments for Column Generation} \label{subsubsec:ret-cg}
We test CG using different pricing heuristics, for solving the LP relaxation of GCP-SC. Specifically, at an iteration of CG, we use a heuristic method to tackle the NP-hard MWISP aiming to quickly find NRC MISs, and add to the RMP at most $n$ NRC MISs in increasing order of their reduced costs. When the heuristic method fails to find any NRC MIS, an exact method, TSM~\citep{jiang2018two}, is used to optimally solve the pricing problem. If the optimal MIS has a negative reduced cost, it is included in the RMP to start the next iteration. Note that the setup of a pricing heuristic and the column-selection strategy are common practices for CG~\citep{lubbecke2005selected}.

We use the MSSP implementation for MWCP described in Section~\ref{subsec:mwcp}. To collect the training data, we use CG to solve LPs on a set of $10$ graphs and record the optimally solved MWISPs encountered during the execution of CG. The graphs are from the graph coloring benchmarks\footnote{\url{https://sites.google.com/site/graphcoloring/files}}, including: ``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". These graphs are selected under two considerations. Firstly, their pricing problems are easy to solve, i.e., collecting the training data is easy. Secondly, they are of different characteristics so that a ML model can learn common knowledge that is applicable to different unseen graphs.

For the compared pricing methods, we include SSSP~\citep{shen2022enhancing}, ACO~\citep{xu2007improved}, Gurobi~\citep{gurobi2018gurobi} and specialized methods to MWCP. These specialized methods include TSM~\citep{jiang2018two}, LSCC+BMS~\citep{wang2016two}, and Fastwclq~\citep{cai2016fast}. For each method, we set a maximum execution time of $30$ seconds for solving a MWISP instance, to prevent one from spending too much time at a particular CG iteration. Note that the exact methods, Gurobi and TSM, are also subject to the cutoff time and treated as `heuristics'. Furthermore, we classify the compared methods into several classes, and the termination criterion for each class is tuned empirically to maximize the performance of CG~\citep{shen2022enhancing}. In particular, the sampling-based methods (i.e., MSSP, SSSP, ACO, and Fastwclq) construct a set of $50n$ MISs at a CG iteration, where $n$ is the number of vertices in a graph. We set MSSP to perform $50$ iterations with $n$ MISs generated at each iteration.

Our test problem instances are generated using the graph coloring benchmarks. A problem instance is a RMP initialized with $10n$ randomly generated columns (or MISs). We generate $24$ problem instances with different random seeds $[1,\cdots,24]$ on a set of $81$ test graphs, resulting in a total number of $1944$ test problem instances. These graphs are used for instance generation, because their initial RMPs can be solved within a reasonable time. That is $500$ seconds under an overall cutoff time of $1800$ seconds in our case. Note that we use the default LP solver in Gurobi~\citep{gurobi2018gurobi} for solving RMPs in CG. All of our experiments are conducted on a server with $32$ CPUs (AMD EPYC Processor, 2245 MHz) and $128$ GB RAM.

\begin{figure}[h!]
	\begin{tikzpicture}
	\begin{groupplot}[group style = {group size = 2 by 1, horizontal sep = 60pt}, height=0.475\textwidth, width=0.475\textwidth, grid style={line width=.1pt, draw=gray!10},major grid style={line width=.2pt,draw=gray!30}, xlabel = \small Time in Seconds, xtick = {0,600,1200,1800}, ticklabel style = {font=\small}, xmajorgrids=true, ymajorgrids=true,  major tick length=0.05cm, minor tick length=0.0cm, legend style={font=\small, column sep = 1pt, legend columns = 4,draw=none}]
   \nextgroupplot[%
    legend to name=group1,
    ylabel = \small \# Solve Instances,
    y label style={at={(axis description cs:-0.12,.5)},anchor=south},
    ]

	\addplot[color=red, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/solving-curve/mssp-cs0.txt};\addlegendentry{\small CG-MSSP}
	\addplot[color=blue, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/solving-curve/sssp-cs0.txt};\addlegendentry{\small CG-SSSP}
	\addplot[color=purple, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/solving-curve/aco.txt};\addlegendentry{\small CG-ACO}
	\addplot[color=orange, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/solving-curve/fastwclq.txt};\addlegendentry{\small CG-Fastwclq}
	\addplot[color=gray, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/solving-curve/lscc.txt};\addlegendentry{\small CG-LSCC}
	\addplot[color=olivegreen, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/solving-curve/tsm.txt};\addlegendentry{\small CG-TSM}
	\addplot[color=black, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/solving-curve/gurobi.txt};\addlegendentry{\small CG-Gurobi}
   \nextgroupplot[%
    ylabel = \small LP Obective,
    y label style={at={(axis description cs:-0.12,.5)},anchor=south},
    ]

	\addplot[color=red, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/lp-curve/mssp-cs0.txt};
	\addplot[color=blue, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/lp-curve/sssp-cs0.txt};
	\addplot[color=purple, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/lp-curve/aco.txt};
	\addplot[color=orange, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/lp-curve/fastwclq.txt};
	\addplot[color=gray, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/lp-curve/lscc.txt};
	\addplot[color=olivegreen, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/lp-curve/tsm.txt};
	\addplot[color=black, line width=0.45mm] table [x=x, y=y, col sep=comma] {cg/small/lp-curve/gurobi.txt};
    \end{groupplot}
    \node at (group c1r1.north) [anchor=north, yshift=1.5cm, xshift=4cm] {\pgfplotslegendfromname{group1}};
	\end{tikzpicture}
\caption{Results for CG with different pricing methods. \textbf{Left:} the number of solved instances out of $1944$ test instances. \textbf{Right:} the objective values of the RMP (the lower the better), averaged over all problem instances using the geometric mean.}
\label{fig:ret_small}
\end{figure}
\paragraph{Results and analysis:}
Figure~\ref{fig:ret_small} shows the results for CG using different pricing heuristics. In the left sub-figure, we report the total number of solved problem instances within $1800$ seconds. In the right sub-figure, we report the progress of CG with respect to the LP objective value, which takes the unsolved problem instances into account as well. It is interesting to observe the following: 1) CG with the ML-based pricing methods (MSSP or SSSP) can make more substantial progress than CG with other methods. This is because ML-based pricing heuristics can generate \emph{many high-quality} columns efficiently. Both the number of NRC columns and the quality of the columns are crucial to the computational efficiency of CG, as noted by previous studies~\citep{lubbecke2005selected}. In at least one of the aspects, ML-based heuristic methods are superior to other pricing methods~\citep{shen2022enhancing}. For example, the performance of CG with LSCC is not competitive because LSCC cannot find a large number of NRC MISs, though LSCC is shown to be very efficient in finding a single (near-)optimal solution~(Figure~\ref{fig:mwc}). 2) A comparison of the ML-based heuristic methods shows that MSSP further improves the performance of CG compared to SSSP. Specifically, CG with MSSP solves many more problem instances than CG with SSSP at any given runtime, as shown in the left sub-figure. This is because MSSP can find better-quality MISs than SSSP for solving MWISP (or MWCP), as evidenced by the results in Figure~\ref{fig:mwc}. However, the results in the right sub-figure of Figure~\ref{fig:ret_small} show that CG-MSSP is marginally better than CG-SSSP in terms of the mean LP objective values. These observations reveal an interesting phenomenon in CG called the tailing-off effect, i.e., the difficulty of CG in converging to the optimal LP solution given near-optimal ones.

\subsection{Experiments for Branch-and-Price}

We test B\&P using CG with different pricing heuristics, for solving the integer problem GCP. We use the B\&P code from an academic mixed-integer-programming solver, SCIP~\citep{gamrath2020scip}. In particular, for solving the pricing problems, the B\&P code uses an efficient greedy search~\citep{mehrotra1996column} as the pricing heuristic and an exact method called $t$-clique. Besides that, this B\&P implementation incorporates some specialized techniques for solving GCPs from previous studies~\citep{mehrotra1996column,malaguti2011exact}, 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. These components form a reasonable B\&P implementation for testing different pricing heuristics. Moreover, we note that additional techniques may further improve the performance of B\&P, such as primal heuristics specialized to CG~\citep{JoncourMSSV10} or methods for stabilizing dual solutions of RMPs~\citep{du1999stabilized}.

We refer to the default setup of B\&P as B\&P-def, and compare it with the B\&P variants that replace the greedy search to one of our ML-based heuristic pricing methods (i.e., MSSP or SSSP), respectively referred to as B\&P-MSSP and B\&P-SSSP. B\&P-SSSP uses a sample size of $10n$ by default~\citep{shen2022enhancing}. To ensure a fair comparison between the ML-based pricing methods, we use B\&P-MSSP to construct the same number of samples evenly divided into $10$ MSSP iterations. From newly generated NRC MISs, we add at most $\theta$ MISs into the RMP in the increasing order of their reduced costs, $\theta=n$ for the root node and $\theta=0.1n$ for child nodes. This is because typically a child node already contain a large number of useful columns inherent from its parent. Hence including a large number of columns is rarely helpful but increases the computational overhead of solving subsequent RMPs.

For each method, we perform $1584$ runs on a set of~$66$ graphs in the graph coloring benchmarks. The excluded graphs are either too easy (all methods can solve them within $10$ seconds) or too hard to solve (all 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 duality gap, defined as
\begin{equation}
\text{duality gap} \coloneqq 100\% \times \frac{upper\_bound - global\_lower\_bound}{upper\_bound},
\end{equation}

\noindent where the upper bound is the objective value of the best found solution and the global lower bound is the smallest lower bound amongst the remaining open tree nodes.

\begin{table}[h!]
    \caption{Results for B\&P with different pricing methods. For each graph, a method solves $24$ problem instances (generated with different random columns in the initial root RMP) and the average results are reported. The solving time is an average over optimally solved GCPs, and `N/A' is placed if a method does not solve any instances on that graph in 8000 seconds. Similarly, the duality gap for a graph is an average over instances whose root LP has been solved to optimality, and `N/A' is placed if the root LP cannot be solved in any run.}
    \label{tab:bp}
    \centering
    \resizebox{\linewidth}{!}{\begin{tabular}{@{}clrr|rrr|rrr|rrr|rrr@{}}
        \toprule
       \multirow{2}{*}{Group} & \multirow{2}{*}{Graph} & \multirow{2}{*}{\# nodes} & \multirow{2}{*}{Density} & \multicolumn{3}{c}{\# optimally Solved} & \multicolumn{3}{c}{Solving time (sec)} & \multicolumn{3}{c}{\# root LP solved} & \multicolumn{3}{c}{Duality gap (\%)} \\
    & & & & MSSP & SSSP & Greedy & MSSP & SSSP & Greedy & MSSP & SSSP & Greedy & MSSP & SSSP & Greedy\\
    \cmidrule(lr){1-1}\cmidrule{2-4}\cmidrule(lr){5-7}\cmidrule(lr){8-10}\cmidrule(lr){11-13}\cmidrule(lr){14-16}

    \multirow{13}{*}{1}  & DSJR500.5 & 486 & 0.972 & \textbf{24} & 17 & 15 & 786.3 & 1370.9 & 782.4 & 24 & 24 & 24 & \textbf{0.0} & 0.2 & 0.3 \\
     & 1-FullIns\_4 & 38 & 0.364 & \textbf{24} & 18 & 23 & \textbf{411.7} & 1289.2 & 675.0  & 24 & 24 & 24 & \textbf{0.0} & 5.0 & 0.8 \\
     & queen9\_9 & 81 & 0.652 & \textbf{24} & 0 & 2 & \textbf{12.6} & N/A & 2204.6 & 24 & 24 & 24 & \textbf{0.0} & 10.0 & 9.2 \\
     & DSJC125.5 & 125 & 0.502 & \textbf{2} & 0 & 0 & \textbf{4905.3} & N/A & N/A & 24 & 24 & 24 & \textbf{9.9} & 15.4 & 15.8 \\
     &  2-Insertions\_3 & 37 & 0.216 & \textbf{24} & 17 & 9 & 4105.3 & 4571.5 & 3448.5 & 24 & 24 & 24 & \textbf{0.0} & 7.3 & 15.6 \\
     & ash608GPIA & 1215 & 0.021 & 24 & 24 & 0 & \textbf{2840.8} & 4899.6 & N/A & 24 & 24 & 0 & 0.0 & 0.0 & N/A \\
     & will199GPIA & 660 & 0.054 & 24 & 24 & 24 & \textbf{54.8} & 70.4 & 128.9 & 24 & 24 & 24 & 0.0 & 0.0 & 0.0 \\
     & 1-Insertions\_4 & 67 & 0.210 & 0 & 0 & 0 & N/A & N/A & N/A & 24 & 24 & 24 & \textbf{20.0} & 37.5 & 40.0 \\
     & myciel5 & 47 & 0.437 & 0 & 0 & 0 & N/A & N/A & N/A & 24 & 24 & 24 & \textbf{16.7} & 22.9 & 32.6 \\
     & DSJC250.9 & 250 & 0.896 & 0 & 0 & 0 & N/A & N/A & N/A & 24 & 24 & 24 & \textbf{2.1} & 2.3 & 4.1 \\
     & r1000.1c & 709 & 0.969 & 0 & 0 & 0 & N/A & N/A & N/A & 24 & 24 & 24 & \textbf{5.9} & 6.2 & 7.1 \\
     & 1-FullIns\_5 & 78 & 0.277 & 0 & 0 & 0 & N/A & N/A & N/A & 24 & 24 & 24 & \textbf{16.7} & 31.9 & 33.3 \\
     & myciel6 & 95 & 0.338 & 0 & 0 & 0 & N/A & N/A & N/A & 24 & 24 & 24 & \textbf{28.6} & 31.0 & 42.9 \\

    \midrule
    \multirow{6}{*}{2}  & le450\_5c & 450 & 0.194 & 13 & \textbf{21} & 1 & 4263.3 & 4024.2 & 2009.9 & 13 & \textbf{21} & 1 & 0.0 & 0.0 & 0.0 \\
    & le450\_5d & 450 & 0.193 & 17 & \textbf{21} & 10 & 4483.0 & 2905.7 & 2126.4 & 21 & \textbf{23} & 11 & 3.2 & 1.5 & 1.5 \\
    & flat300\_26\_0 & 300 & 0.965 & 24 & 24 & 24 & 361.6 & \textbf{223.0} & 1910.6 & 24 & 24 & 24 & 0.0 & 0.0 & 0.0 \\
    & school1\_nsh & 326 & 0.547 & 24 & 24 & 24 & 685.8 & \textbf{627.0} & 759.5 & 24 & 24 & 24 & 0.0 & 0.0 & 0.0 \\
    & ash331GPIA & 661 & 0.038 & 24 & 24 & 24 & 59.1 & \textbf{53.1} & 307.4 & 24 & 24 & 24 & 0.0 & 0.0 & 0.0 \\
    & queen16\_16 & 256 & 0.387 &  0 & 0  & 0  & N/A & N/A & N/A & 23 & \textbf{24} & 0 & 11.1 & 11.1 & N/A \\

    \midrule
    \multirow{10}{*}{3}  & flat300\_20\_0 & 300 & 0.953 & 24 & 24 & 24 & 84.9 & 84.9 & 343.6 & 24 & 24 & 24 & 0.0 & 0.0 & 0.0 \\
    & school1 & 355 & 0.603 & 24 & 24 & 24 & 51.7 & 52.0 & 1966.9 & 24 & 24 & 24 & 0.0 & 0.0 & 0.0 \\
    & DSJC1000.9 & 1000 & 0.900 &  0 & 0  & 0  & N/A & N/A & N/A & 24 & 24 & 17 & 11.5 & 11.5 & 11.5 \\
    & le450\_25d & 433 & 0.366 &  0 & 0  & 0  & N/A & N/A & N/A & 24 & 24 & 0 & 10.7 & 10.7 & N/A \\
    & le450\_25c & 435 & 0.362 &  0 & 0  & 0  & N/A & N/A & N/A & 24 & 24 & 0 & 10.7 & 10.7 & N/A \\
    & le450\_15b & 410 & 0.187 &  0 & 0  & 0  & N/A & N/A & N/A & 24 & 24 & 0 & 6.3 & 6.3 & N/A \\
    &  qg.order40 & 1600 & 0.098 &  0 & 0  & 0  & N/A & N/A & N/A & 24 & 24 & 0 & 2.4 & 2.4 & N/A \\
    & le450\_15a & 407 & 0.189 &  0 & 0  & 0  & N/A & N/A & N/A & 24 & 24 & 0 & 6.3 & 6.3 & N/A \\
    & wap06a & 703 & 0.288 &  0 & 0  & 0  & N/A & N/A & N/A & 24 & 24 & 0 & 4.8 & 4.8 & N/A \\
    & queen15\_15 & 225 & 0.411 &  0 & 0  & 0  & N/A & N/A & N/A & 24 & 24 & 0 & 11.8 & 11.8 & N/A \\

    \midrule
    \multirow{5}{*}{4} & 3-Insertions\_5 & 1406 & 0.020 &  0 & 0  & 0  & N/A & N/A & N/A & 1 & 0 & \textbf{24} & 50.0 & N/A & 50.0 \\
    & 2-Insertions\_5 & 597 & 0.044 &  0 & 0  & 0  & N/A & N/A & N/A & 1 & 1 & \textbf{24} & 50.0 & 50.0 & 50.0 \\
    & 4-Insertions\_4 & 475 & 0.032 &  0 & 0  & 0  & N/A & N/A & N/A & 23 & 15 & \textbf{24} & 40.0 & 40.0 & 40.0 \\
    & r1000.5 & 966 & 0.477 &  0 & 0  & 0  & N/A & N/A & N/A & 24 & 24 & 24 & 9.3 & 7.8 & \textbf{1.3} \\
    & 1-Insertions\_5 & 202 & 0.121 &  0 & 0  & 0  & N/A & N/A & N/A & 24 & 24 & 24 & 34.0 & 43.8 & \textbf{33.3} \\
        \bottomrule
    \end{tabular}}
\end{table}

\paragraph{Results and analysis:}
Table~\ref{tab:bp} shows the results for B\&P. According to the performances of B\&P using different pricing heuristics, the test graphs are divided into four groups and are described respectively. Group~1 includes $13$ graphs where B\&P with MSSP outperforms others, and the improvements are substantial on several graphs. Specifically, B\&P-MSSP solves GCP on `queen9\_9' in all runs with an average time of $12.6$ seconds and on `DSJC125.5' in two runs, whereas other methods cannot solve any with a cutoff time of $8000$ seconds; B\&P-MSSP reduces the duality gap almost by half compared to other methods for the graphs `1-Insertions\_4' and `1-FullIns-5'. Notably, B\&P-MSSP also significantly improves on B\&P-SSSP. This is expected as B\&P repeatedly invokes CG at each node of the B\&P process, and CG-MSSP is more efficient than CG-SSSP in general (see Section~\ref{subsubsec:ret-cg}). Group 2 consists of $6$ graphs on which B\&P with SSSP outperforms others. In particular, B\&P-SSSP is considerably better than B\&P-MSSP on the two graphs with the prefix `le450-5'. For these graphs, the ML predictions by MSSP may have converged fast to local optima, resulting in a less diverse set of sampled MISs. Group~3 shows the results on $10$ graphs, where the performances of B\&P with the two ML-based pricing heuristics are comparable and are much better than B\&P with its default greedy search. This is because for these graphs CG can be sped up substantially by including multiple MISs at each iteration. Both ML-based pricing methods can generate many high-quality MISs, whereas the default greedy heuristic can only generate a single MIS per CG iteration. In contrast, results in Group~4 show that B\&P with the ML-based pricing methods do not outperform B\&P-def. This is because, for these graphs, the exact method can efficiently solve the pricing problems already, and competitive performance of CG can be obtained by disabling the use of a heuristic method~\citep{shen2022enhancing}. Lastly, we note that on some test graphs the results for the three methods are comparable, and further experimental results can be found in the Supplementary Materials.

\section{Discussions}

\section{Conclusion}
\label{sec:conclusion}
This paper proposes a multi-shot solution prediction (MSSP) framework that enables a ML model to predict optimal solutions in multiple shots. We instantiate MSSP using support vector machine with linear kernel (linear-SVM) and a probabilistic sampling (PS) method, where the linear-SVM produces solution predictions to guide the PS method for finding high-quality solutions, and those solutions are fed back to linear-SVM via statistical features to further improve the quality of the prediction.

The efficacy of MSSP is demonstrated in two scenarios. Firstly, on three NP-hard combinatorial problems, our experiments show that the proposed MSSP can significantly improve the prediction quality of linear-SVM, and the more accurate ML prediction can in turn guide the search to find better-quality solutions efficiently. Moreover, MSSP can substantially outperform a range of methods, such as ML-based approaches, ant colony optimization, and a mixed-integer-programming solver. Secondly, on the graph coloring problem (GCP), our experiments show that using MSSP as a pricing heuristic can boost column generation (CG) for solving the linear-programming relaxation of GCP. As a result, the performance of the branch-and-price (B\&P) method can be improved substantially.

In future work, we would like to test CG using our MSSP as pricing heuristics on more combinatorial optimization problems such as vehicle routing problems. Our goal is to develop a \emph{generic} ML-based pricing heuristic to speed up CG and B\&P for solving the Dantzig-Wolfe reformulation of combinatorial optimization problems.

redpajama_set_name: RedPajamaArXiv