# Exact Combinatorial Optimization
with Graph Convolutional Neural Networks

Maxime Gasse
Mila, Polytechnique Montréal
<EMAIL_ADDRESS>
&Didier Chételat
Polytechnique Montréal
<EMAIL_ADDRESS>
&Nicola Ferroni
University of Bologna
<EMAIL_ADDRESS>
&Laurent Charlin
Mila, HEC Montréal
<EMAIL_ADDRESS>
&Andrea Lodi
Mila, Polytechnique Montréal
<EMAIL_ADDRESS>

###### Abstract

Combinatorial optimization problems are typically tackled by the branch-and-
bound paradigm. We propose a new graph convolutional neural network model for
learning branch-and-bound variable selection policies, which leverages the
natural variable-constraint bipartite graph representation of mixed-integer
linear programs. We train our model via imitation learning from the strong
branching expert rule, and demonstrate on a series of hard problems that our
approach produces policies that improve upon state-of-the-art machine-learning
methods for branching and generalize to instances significantly larger than
seen during training. Moreover, we improve for the first time over expert-
designed branching rules implemented in a state-of-the-art solver on large
problems. Code for reproducing all the experiments can be found at
https://github.com/ds4dm/learn2branch.

##  1 Introduction

Combinatorial optimization aims to find optimal configurations in discrete
spaces where exhaustive enumeration is intractable. It has applications in
fields as diverse as electronics, transportation, management, retail, and
manufacturing [42], but also in machine learning, such as in structured
prediction and maximum a posteriori inference [51, 34, 49]. Such problems can
be extremely difficult to solve, and in fact most classical NP-hard computer
science problems are examples of combinatorial optimization. Nonetheless,
there exists a broad range of exact combinatorial optimization algorithms,
which are guaranteed to find an optimal solution despite a worst-case
exponential time complexity [52]. An important property of such algorithms is
that, when interrupted before termination, they can usually provide an
intermediate solution along with an optimality bound, which can be a valuable
information in theory and in practice. For example, after one hour of
computation, an exact algorithm may give the guarantee that the best solution
found so far lies within 2% of the optimum, even without knowing what the
actual optimum is. This quality makes exact methods appealing and practical,
and as such they constitute the core of modern commercial solvers.

In practice, most combinatorial optimization problems can be formulated as
mixed-integer linear programs (MILPs), in which case branch-and-bound (B&B)
[35] is the exact method of choice. Branch-and-bound recursively partitions
the solution space into a search tree, and computes relaxation bounds along
the way to prune subtrees that provably cannot contain an optimal solution.
This iterative process requires sequential decision-making, such as _node
selection_ : selecting the next node to evaluate, and _variable selection_ :
selecting the variable by which to partition the node’s search space [41].
This decision process traditionally follows a series of hard-coded heuristics,
carefully designed by experts to minimize the average solving time on a
representative set of MILP instances [21]. However, in many contexts it is
common to repeatedly solve similar combinatorial optimization problems, e.g.,
day-to-day production planning and lot-sizing problems [44], which may
significantly differ from the set of instances on which B&B algorithms are
typically evaluated. It is then appealing to use statistical learning for
tuning B&B algorithms automatically for a desired class of problems. However,
this line of work raises two challenges. First, it is not obvious how to
encode the state of a MILP B&B decision process [4], especially since both
search trees and integer linear programs can have a variable structure and
size. Second, it is not clear how to formulate a model architecture that leads
to rules which can generalize, at least to similar instances but also ideally
to instances larger than seen during training.

In this work we propose to address the above challenges by using graph
convolutional neural networks. More precisely, we focus on variable selection,
also known as the _branching problem_ , which lies at the core of the B&B
paradigm yet is still not well theoretically understood [41], and adopt an
imitation learning strategy to learn a fast approximation of _strong
branching_ , a high-quality but expensive branching rule. While such an idea
is not new [30, 4, 24], we propose to address the learning problem in a novel
way, through two contributions. First, we propose to encode the branching
policies into a graph convolutional neural network (GCNN), which allows us to
exploit the natural bipartite graph representation of MILP problems, thereby
reducing the amount of manual feature engineering. Second, we approximate
strong branching decisions by using behavioral cloning with a cross-entropy
loss, a less difficult task than predicting strong branching scores [4] or
rankings [30, 24]. We evaluate our approach on four classes of NP-hard
problems, namely set covering, combinatorial auction, capacitated facility
location and maximum independent set. We compare against the previously
proposed approaches of Khalil et al. [30], Alvarez et al. [4] and Hansknecht
et al. [24], as well as against the default hybrid branching rule in SCIP
[20], a modern open-source solver. The results show that our choice of model,
state encoding, and training procedure leads to policies that can offer a
substantial improvement over traditional branching rules, and generalize well
to larger instances than those used in training.

In Section 2, we review the broader literature of works that use statistical
learning for branching. In Section 3, we formally introduce the B&B framework,
and formulate the branching problem as a Markov decision process. In Section
4, we present our state representation, model, and training procedure for
addressing the branching problem. Finally, we discuss experimental results in
Section 5.

##  2 Related work

First steps towards statistical learning of branching rules in B&B were taken
by Khalil et al. [30], who learn a branching rule customized to a single
instance during the B&B process, as well as Alvarez et al. [4] and Hansknecht
et al. [24] who learn a branching rule offline on a collection of similar
instances, in a fashion similar to us. In each case a branching policy is
learned by imitation of the strong branching expert, although with a
differently formulated learning problem. Namely, Khalil et al. [30] and
Hansknecht et al. [24] treat it as a ranking problem and learn a partial
ordering of the candidates produced by the expert, while Alvarez et al. [4]
treat it as a regression problem and learn directly the strong branching
scores of the candidates. In contrast, we treat it as a classification problem
and simply learn from the expert decisions, which allows imitation from
experts that don’t rely on branching scores or orderings. These works also
differ from ours in three other key aspects. First, they rely on extensive
feature engineering, which is reduced by our graph convolutional neural
network approach. Second, they do not evaluate generalization ability to
instances larger than seen during training, which we propose to do. Finally,
in each case performance was evaluated on a simplified solver, whereas we
compare, for the first time and favorably, against a full-fledged solver with
primal heuristics, cuts and presolving activated. We compare against these
approaches in Section 5.

Other works have considered using graph convolutional neural networks in the
context of approximate combinatorial optimization, where the objective is to
find good solutions quickly, without seeking any optimality guarantees. The
first work of this nature was by Khalil et al. [31], who proposed a GCNN model
for learning greedy heuristics on several collections of combinatorial
optimization problems defined on graphs. This was followed by Selsam et al.
[47], who proposed a recurrent GCNN model, NeuroSAT, which can be interpreted
as an approximate SAT solver when trained to predict satisfiability. Such
works provide additional evidence that GCNNs can effectively capture
structural characteristics of combinatorial optimization problems.

Other works consider using machine learning to improve variable selection in
branch-and-bound, without directly learning a branching policy. Di Liberto et
al. [15] learn a clustering-based classifier to pick a variable selection rule
at every branching decision up to a certain depth, while Balcan et al. [8] use
the fact that many variable selection rules in B&B explicitly score the
candidate variables, and propose to learn a weighting of different existing
scores to combine their strengths. Other works learn variable selection
policies, but for algorithms less general than B&B. Liang et al. [39] learn a
variable selection policy for SAT solvers using a bandit approach, and
Lederman et al. [36] extend their work by taking a reinforcement learning
approach with graph convolutional neural networks. Unlike our approach, these
works are restricted to conflict-driven clause learning methods in SAT
solvers, and cannot be readily extended to B&B methods for arbitrary mixed-
integer linear programs. In the same vein, Balunovic et al. [9] learn by
imitation learning a variable selection procedure for SMT solvers that
exploits specific aspects of this type of solver.

Finally, researchers have also focused on learning other aspects of B&B
algorithms than variable selection. He et al. [25] learn a node selection
heuristic by imitation learning of the oracle procedure that expands the node
whose feasible set contains the optimal solution, while Song et al. [48] learn
node selection and pruning heuristics by imitation learning of shortest paths
to good feasible solutions, and Khalil et al. [32] learn primal heuristics for
B&B algorithms. Those approaches are complementary with our work, and could in
principle be combined to further improve solver performance. More generally,
many authors have proposed machine learning approaches to fine-tune exact
optimization algorithms, not necessarily for MILPs in general. A recent survey
is provided by Bengio et al. [10].

##  3 Background

###  3.1 Problem definition

A mixed-integer linear program is an optimization problem of the form

$\displaystyle\operatorname*{arg\,min}_{\mathbf{x}}\left\\{\mathbf{c}^{\top}\mathbf{x}\mid\mathbf{A}\mathbf{x}\leq\mathbf{b},\enskip\mathbf{l}\leq\mathbf{x}\leq\mathbf{u},\enskip\mathbf{x}\in\mathbb{Z}^{p}\times\mathbb{R}^{n-p}\right\\}\text{,}$
(1)

where $\mathbf{c}\in\mathbb{R}^{n}$ is called the objective coefficient
vector, $\mathbf{A}\in\mathbb{R}^{m\times n}$ the constraint coefficient
matrix, $\mathbf{b}\in\mathbb{R}^{m}$ the constraint right-hand-side vector,
$\mathbf{l},\mathbf{u}\in\mathbb{R}^{n}$ respectively the lower and upper
variable bound vectors, and $p\leq n$ the number of integer variables. Under
this representation, the size of a MILP is typically measured by the number of
rows ($m$) and columns ($n$) of the constraint matrix. By relaxing the
integrality constraint, one obtains a continuous linear program (LP) whose
solution provides a lower bound to (1), and can be solved efficiently using,
for example, the simplex algorithm. If a solution to the LP relaxation
respects the original integrality constraint, then it is also a solution to
(1). If not, then one may decompose the LP relaxation into two sub-problems,
by splitting the feasible region according to a variable that does not respect
integrality in the current LP solution $\mathbf{x}^{\star}$,

$\displaystyle x_{i}\leq\lfloor x_{i}^{\star}\rfloor\lor x_{i}\geq\lceil
x_{i}^{\star}\rceil,\quad\exists i\leq p\mid
x_{i}^{\star}\not\in\mathbb{Z}\text{,}$ (2)

where $\lfloor.\rfloor$ and $\lceil.\rceil$ respectively denote the floor and
ceil functions. In practice, the two sub-problems will only differ from the
parent LP in the variable bounds for $x_{i}$, which get updated to
$u_{i}=\lfloor x_{i}^{\star}\rfloor$ in the left child and $l_{i}=\lceil
x_{i}^{\star}\rceil$ in the right child.

The branch-and-bound algorithm [52, Ch. II.4], in its simplest formulation,
repeatedly performs this binary decomposition, giving rise to a search tree.
By design, the best LP solution in the leaf nodes of the tree provides a lower
bound to the original MILP, whereas the best integral LP solution (if any)
provides an upper bound. The solving process stops whenever both the upper and
lower bounds are equal or when the feasible regions do not decompose anymore,
thereby providing a certificate of optimality or infeasibility, respectively.

###  3.2 Branching rules

A key step in the B&B algorithm is selecting a fractional variable to branch
on in (2), which can have a very significant impact on the size of the
resulting search tree [2]. As such, branching rules are at the core of modern
combinatorial optimization solvers, and have been the focus of extensive
research [40, 43, 1, 17]. So far, the branching strategy consistently
resulting in the smallest B&B trees is _strong branching_ [5]. It does so by
computing the expected bound improvement for each candidate variable before
branching, which unfortunately requires the solution of two LPs for every
candidate. In practice, running strong branching at every node is prohibitive,
and modern B&B solvers instead rely on _hybrid branching_ [3, 1] which
computes strong branching scores only at the beginning of the solving process
and gradually switches to simpler heuristics such as: the conflict score (in
the original article), the pseudo-cost [43] or a hand-crafted combination of
the two. For a more extensive discussion of B&B branching strategies in MILP,
the reader is referred to Achterberg et al. [3].

###  3.3 Markov decision process formulation

As remarked by He et al. [25], the sequential decisions made during B&B can be
assimilated to a Markov decision process [26]. Consider the solver to be the
environment, and the brancher the agent. At the $t^{\text{th}}$ decision the
solver is in a state $\mathbf{s}_{t}$, which comprises the B&B tree with all
past branching decisions, the best integer solution found so far, the LP
solution of each node, the currently focused leaf node, as well as any other
solver statistics (such as, for example, the number of times every primal
heuristic has been called). The brancher then selects a variable
$\mathbf{a}_{t}$ among all fractional variables
$\mathcal{A}(\mathbf{s}_{t})\subseteq\\{1,\dots,p\\}$ at the currently focused
node, according to a policy $\pi(\mathbf{a}_{t}\,|\,\mathbf{s}_{t})$. The
solver in turn extends the B&B tree, solves the two child LP relaxations, runs
any internal heuristic, prunes the tree if warranted, and finally selects the
next leaf node to split. We are then in a new state $\mathbf{s}_{t+1}$, and
the brancher is called again to take the next branching decision. This
process, illustrated in Figure 1, continues until the instance is solved,
i.e., until there are no leaf node left for branching.

As a Markov decision process, B&B is episodic, where each episode amounts to
solving a MILP instance. Initial states correspond to an instance being
sampled among a group of interest, while final states mark the end of the
optimization process. The probability of a trajectory
$\tau=(\mathbf{s}_{0},\dots,\mathbf{s}_{T})\in\mathcal{T}$ then depends on
both the branching policy $\pi$ and the remaining components of the solver,

$p_{\pi}(\tau)=p(\mathbf{s}_{0})\prod_{t=0}^{T-1}\sum_{\mathbf{a}\in\mathcal{A}(\mathbf{s}_{t})}\pi(\mathbf{a}\,|\,\mathbf{s}_{t})p(\mathbf{s}_{t+1}\,|\,\mathbf{s}_{t},\mathbf{a})\text{.}$

A natural approach to find good branching policies is reinforcement learning,
with a carefully designed reward function. However, this raises several key
issues which we circumvent by adopting an imitation learning scheme, as
discussed next.

$\mathbf{s}_{t}$$x_{7}\leq 0$$x_{7}\geq 1$$x_{1}\leq 2$$x_{1}\geq
3$$\mathcal{A}(\mathbf{s}_{t})=\\{1,3,4\\}$$\mathbf{a}_{t}=4$$\mathbf{s}_{t+1}$$x_{7}\leq
0$$x_{7}\geq 1$$x_{1}\leq 2$$x_{1}\geq 3$$x_{4}\leq-2$$x_{4}\geq-1$

Figure 1: B&B variable selection as a Markov decision process. On the left, a
state $\mathbf{s}_{t}$ comprised of the branch-and-bound tree, with a leaf
node chosen by the solver to be expanded next (in pink). On the right, a new
state $\mathbf{s}_{t+1}$ resulting from branching on the variable
$\mathbf{a}_{t}=x_{4}$.

##  4 Methodology

We now describe our approach for tackling the B&B variable selection problem
in MILPs, where we use imitation learning and a dedicated graph convolutional
neural network model. As the B&B variable selection problem can be formulated
as a Markov decision process, a natural way of training a policy would be
reinforcement learning [50]. However, this approach runs into many issues.
Notably, as episode length is proportional to performance, and randomly
initialized policies perform poorly, standard reinforcement learning
algorithms are usually so slow early in training as to make total training
time prohibitively long. Moreover, once the initial state corresponding to an
instance is selected, the rest of the process is instance-specific, and so the
Markov decision processes tend to be extremely large. In this work we choose
instead to learn directly from an expert branching rule, an approach usually
referred to as imitation learning [27].

###  4.1 Imitation learning

We train by behavioral cloning [45] using the strong branching rule, which
suffers a high computational cost but usually produces the smallest B&B trees,
as mentioned in Section 3.2. We first run the expert on a collection of
training instances of interest, record a dataset of expert state-action pairs
$\mathcal{D}=\\{(\mathbf{s}_{i},\mathbf{a}^{\star}_{i})\\}_{i=1}^{N}$, and
then learn our policy by minimizing the cross-entropy loss

$\mathcal{L}(\theta)=-\frac{1}{N}\sum_{(\mathbf{s},\mathbf{a}^{*})\in\mathcal{D}}\log\pi_{\theta}(\mathbf{a}^{*}\,|\,\mathbf{s}).$
(3)

$\mathbf{v}_{1}$$\mathbf{v}_{2}$$\mathbf{v}_{3}$$\mathbf{c}_{1}$$\mathbf{c}_{2}$$\mathbf{e}_{1,1}$$\mathbf{e}_{1,3}$$\mathbf{e}_{1,2}$$\mathbf{e}_{2,3}$

$\mathbf{C}$$\mathbf{E}$$\mathbf{V}$$\mathbf{C}^{1}$$\mathbf{V}^{1}$$\mathbf{C}^{2}$$\mathbf{V}^{2}$$\pi(\mathbf{x})$$n\times
d$$m\times n\times e$$m\times c$$n\times 64$$m\times 64$$m\times 64$$n\times
64$$n\times
1$initialembedding$\mathcal{C}$-sideconvolution$\mathcal{V}$-sideconvolutionfinalembedding\+
softmax

Figure 2: Left: our bipartite state representation
$\mathbf{s}_{t}=(\mathcal{G},\mathbf{C},\mathbf{E},\mathbf{V})$ with $n=3$
variables and $m=2$ constraints. Right: our bipartite GCNN architecture for
parametrizing our policy $\pi_{\theta}(\mathbf{a}\mid\mathbf{s}_{t})$.

###  4.2 State encoding

We encode the state $\mathbf{s}_{t}$ of the B&B process at time $t$ as a
bipartite graph with node and edge features
$(\mathcal{G},\mathbf{C},\mathbf{E},\mathbf{V})$, described in Figure 2
(Left). On one side of the graph are nodes corresponding to the constraints in
the MILP, one per row in the current node’s LP relaxation, with
$\mathbf{C}\in\mathbb{R}^{m\times c}$ their feature matrix. On the other side
are nodes corresponding to the variables in the MILP, one per LP column, with
$\mathbf{V}\in\mathbb{R}^{n\times d}$ their feature matrix. An edge
$(i,j)\in\mathcal{E}$ connects a constraint node $i$ and a variable node $j$
if the latter is involved in the former, that is if $\mathbf{A}_{ij}\not=0$,
and $\mathbf{E}\in\mathbb{R}^{m\times n\times e}$ represents the (sparse)
tensor of edge features. Note that under mere restrictions in the B&B solver
(namely, by enabling cuts only at the root node), the graph structure is the
same for all LPs in the B&B tree, which reduces the cost of feature
extraction. The exact features attached to the graph are described in the
supplementary materials. We note that this is really only a subset of the
solver state, which technically turns the process into a partially-observable
Markov decision process [6], but also that excellent variable selection
policies such as strong branching are able to do well despite relying only on
a subset of the solver state as well.

###  4.3 Policy parametrization

We parametrize our variable selection policy
$\pi_{\theta}(\mathbf{a}|\mathbf{s}_{t})$ as a graph convolutional neural
network [23, 46, 12]. Such models, also known as message-passing neural
networks [19], are extensions of convolutional neural networks from grid-
structured data (as in images or sounds) to arbitrary graphs. They have been
successfully applied to a variety of machine learning tasks with graph-
structured inputs, such as prediction of molecular properties [16, 19],
program verification [38], and document classification in citation networks
[33]. Graph convolutions exhibit many properties which make them a natural
choice for graph-structured data in general, and MILP problems in particular:
1) they are well-defined no matter the input graph size; 2) their
computational complexity is directly related to the density of the graph,
which makes it an ideal choice for processing typically sparse MILP problems;
and 3) they are permutation-invariant, that is they will always produce the
same output no matter the order in which the nodes are presented.

Our model takes as input our bipartite state representation
$\mathbf{s}_{t}=(\mathcal{G},\mathbf{C},\mathbf{V},\mathbf{E})$ and performs a
single graph convolution, in the form of two interleaved half-convolutions. In
detail, because of the bipartite structure of the input graph, our graph
convolution can be broken down into two successive passes, one from variable
to constraints and one from constraints to variables. These passes take the
form

$\mathbf{c}_{i}\leftarrow\mathbf{f}_{\mathcal{C}}\Big{(}\mathbf{c}_{i},\sum_{j}^{(i,j)\in\mathcal{E}}\mathbf{g}_{\mathcal{C}}\left(\mathbf{c}_{i},\mathbf{v}_{j},\mathbf{e}_{i,j}\right)\Big{)}\text{,}\qquad\mathbf{v}_{j}\leftarrow\mathbf{f}_{\mathcal{V}}\Big{(}\mathbf{v}_{j},\sum_{i}^{(i,j)\in\mathcal{E}}\mathbf{g}_{\mathcal{V}}\left(\mathbf{c}_{i},\mathbf{v}_{j},\mathbf{e}_{i,j}\right)\Big{)}$
(4)

for all $i\in\mathcal{C}$, $j\in\mathcal{V}$, where
$\mathbf{f}_{\mathcal{C}}$, $\mathbf{f}_{\mathcal{V}}$,
$\mathbf{g}_{\mathcal{C}}$ and $\mathbf{g}_{\mathcal{V}}$ are 2-layer
perceptrons with relu activation functions. Following this graph-convolution
layer, we obtain a bipartite graph with the same topology as the input, but
with potentially different node features, so that each node now contains
information from its neighbors. We obtain our policy by discarding the
constraint nodes and applying a final 2-layer perceptron on variable nodes,
combined with a masked softmax activation to produce a probability
distribution over the candidate branching variables (i.e., the non-fixed LP
variables). The right side of Figure 2 provides an overview of our
architecture.

#### Prenorm layers

In the literature of GCNN, it is common to normalize each convolution
operation by the number of neighbours [33]. As noted by Xu et al. [53] this
might result in a loss of expressiveness, as the model then becomes unable to
perform a simple counting operation (e.g., in how many constraints does a
variable appears). Therefore we opt for un-normalized convolutions. However,
this introduces a weight initialization issue. Indeed, weight initialization
in standard CNNs relies on the number of input units to normalize the initial
weights [22], which in a GCNN is unknown beforehand and depends on the
dataset. To overcome this issue and stabilize the learning procedure, we adopt
a simple affine transformation
$\mathbf{x}\leftarrow(\mathbf{x}-\beta)/\sigma$, which we call a _prenorm_
layer, applied right after the summation in (4). The $\beta$ and $\sigma$
parameters are initialized with respectively the empirical mean and standard
deviation of $\mathbf{x}$ on the training dataset, and fixed once and for all
before the actual training. Adopting both un-normalized convolutions and this
pre-training procedure improves our generalization performance on larger
problems, as will be shown in Section 5.

##  5 Experiments

We now present a comparative experiment against three competing machine
learning approaches and SCIP’s default branching rule to assess the value of
our approach, as well as an ablation study to validate our architectural
choices. Code for reproducing all the experiments can be found at
https://github.com/ds4dm/learn2branch.

###  5.1 Setup

#### Benchmarks

We evaluate our approach on four NP-hard problem benchmarks. Our first
benchmark is comprised of set covering instances generated following Balas and
Ho [7], with 1,000 columns. We train and test on instances with 500 rows, and
we evaluate on instances with 500 (Easy), 1,000 (Medium) and 2,000 (Hard)
rows. Our second benchmark is comprised of combinatorial auction instances,
generated following the arbitrary relationships procedure of Leyton-Brown et
al. [37, Section 4.3]. We train and test on instances with 100 items for 500
bids, and we evaluate on instances with 100 items for 500 bids (Easy), 200
items for 1,000 bids (Medium) and 300 items for 1,500 bids (Hard). Our third
benchmark is comprised of capacitated facility location instances generated
following Cornuejols et al. [14], with 100 facilities. We train and test on
instances with 100 customers, and we evaluate on instances with 100 (Easy),
200 (Medium) and 400 (Hard) customers. Finally, our fourth benchmark is
comprised of maximum independent set instances on Erdős-Rényi random graphs,
generated following the procedure of Bergman et al. [11, 4.6.4] with affinity
set to $4$. We train and test on instances with graphs of 500 nodes, and we
evaluate on instances with 500 (Easy), 1000 (Medium) and 1500 nodes (Hard).
These four benchmarks were chosen because they are challenging for state-of-
the-art solvers, but also representative of the types of integer programming
problems encountered in practice. In particular, set covering problems capture
the quintessence of integer linear programming, since column generation
formulations can be written for virtually any difficult discrete optimization
problem. Throughout all experiments we use SCIP 6.0.1 as the backend solver,
with a time limit of 1 hour. Following Karzan et al. [29], Fischetti and
Monaci [17] and Khalil et al. [30] we allow cutting plane generation at the
root node only, and deactivate solver restarts. All other SCIP parameters are
kept to default so as to make comparisons as fair and reproducible as
possible.

#### Baselines

We compare against a human-designed state-of-the-art branching rule:
reliability pseudocost (rpb), a variant of hybrid branching [1] which is used
by default in SCIP. For completeness, we report as well the performance of
full strong branching (fsb), our slow expert. We also compare against three
machine learning branchers: the learning-to-score approach of Alvarez et al.
[4] (trees) based on an ExtraTrees [18] model, as well as the learning-to-rank
approaches from Khalil et al. [30] (svmrank) and Hansknecht et al. [24]
(lmart), based on an SVMrank [28] and a LambdaMART [13] model, respectively.
The trees model uses variable-wise features from our bipartite state, obtained
by concatenating the variable’s node features with edge and constraint node
features statistics over its neighborhood. The svmrank and lmart models both
use the original features proposed by Khalil et al. [30], which we re-
implemented within SCIP. More training details for each machine learning
method can be found in the supplementary materials.

Table 1: Imitation learning accuracy on the test sets. | Set Covering | Combinatorial Auction | Capacitated Facility Location | Maximum Independent Set
---|---|---|---|---
model | acc@1 | acc@5 | acc@10 | acc@1 | acc@5 | acc@10 | acc@1 | acc@5 | acc@10 | acc@1 | acc@5 | acc@10
trees |  $51.8$$\pm$ | $0.3$ |  $80.5$$\pm$ | $0.1$ |  $91.4$$\pm$ | $0.2$ |  $52.9$$\pm$ | $0.3$ |  $84.3$$\pm$ | $0.1$ |  $94.1$$\pm$ | $0.1$ |  $63.0$$\pm$ | $0.4$ |  $97.3$$\pm$ | $0.1$ |  $99.9$$\pm$ | $0.0$ |  $30.9$$\pm$ | $0.4$ |  $47.4$$\pm$ | $0.3$ |  $54.6$$\pm$ | $0.3$
svmrank |  $57.6$$\pm$ | $0.2$ |  $84.7$$\pm$ | $0.1$ |  $94.0$$\pm$ | $0.1$ |  $57.2$$\pm$ | $0.2$ |  $86.9$$\pm$ | $0.2$ |  $95.4$$\pm$ | $0.1$ |  $67.8$$\pm$ | $0.1$ |  $98.1$$\pm$ | $0.1$ |  $99.9$$\pm$ | $0.0$ |  $48.0$$\pm$ | $0.6$ |  $69.3$$\pm$ | $0.2$ |  $78.1$$\pm$ | $0.2$
lmart |  $57.4$$\pm$ | $0.2$ |  $84.5$$\pm$ | $0.1$ |  $93.8$$\pm$ | $0.1$ |  $57.3$$\pm$ | $0.3$ |  $86.9$$\pm$ | $0.2$ |  $95.3$$\pm$ | $0.1$ |  $68.0$$\pm$ | $0.2$ |  $98.0$$\pm$ | $0.0$ |  $99.9$$\pm$ | $0.0$ |  $48.9$$\pm$ | $0.3$ |  $68.9$$\pm$ | $0.4$ |  $77.0$$\pm$ | $0.5$
gcnn |  $65.5$$\pm$ | $0.1$ |  $92.4$$\pm$ | $0.1$ |  $98.2$$\pm$ | $0.0$ |  $61.6$$\pm$ | $0.1$ |  $91.0$$\pm$ | $0.1$ |  $97.8$$\pm$ | $0.1$ |  $71.2$$\pm$ | $0.2$ |  $98.6$$\pm$ | $0.1$ |  $99.9$$\pm$ | $0.0$ |  $56.5$$\pm$ | $0.2$ |  $80.8$$\pm$ | $0.3$ |  $89.0$$\pm$ | $0.1$
Table 2: Policy evaluation on separate instances in terms of solving time, number of wins (fastest method) over number of solved instances, and number of resulting B&B nodes (lower is better). For each problem, the models are trained on easy instances only. See Section 5.1 for definitions. | Easy | Medium | Hard
---|---|---|---
Model | Time | Wins | Nodes | Time | Wins | Nodes | Time | Wins | Nodes
fsb |  $17.30$ $\pm$ |  $6.1$% |  $0$ / | $100$ |  $17$ $\pm$ |  $13.7$% |  $411.34$ $\pm$ |  $4.3$% |  $0$ / | $90$ |  $171$ $\pm$ |  $6.4$% |  $3600.00$ $\pm$ |  $0.0$% |  $0$ / | $0$ |  n/a $\pm$ | n/a%
rpb |  $8.98$ $\pm$ |  $4.8$% |  $0$ / | $100$ |  $54$ $\pm$ |  $20.8$% |  $60.07$ $\pm$ |  $3.7$% |  $0$ / | $100$ |  $1741$ $\pm$ |  $7.9$% |  $1677.02$ $\pm$ |  $3.0$% |  $4$ / | $65$ |  $47\,299$ $\pm$ |  $4.9$%
trees |  $9.28$ $\pm$ |  $4.9$% |  $0$ / | $100$ |  $187$ $\pm$ |  $9.4$% |  $92.47$ $\pm$ |  $5.9$% |  $0$ / | $100$ |  $2187$ $\pm$ |  $7.9$% |  $2869.21$ $\pm$ |  $3.2$% |  $0$ / | $35$ |  $59\,013$ $\pm$ |  $9.3$%
svmrank |  $8.10$ $\pm$ |  $3.8$% |  $1$ / | $100$ |  $165$ $\pm$ |  $8.2$% |  $73.58$ $\pm$ |  $3.1$% |  $0$ / | $100$ |  $1915$ $\pm$ |  $3.8$% |  $2389.92$ $\pm$ |  $2.3$% |  $0$ / | $47$ |  $42\,120$ $\pm$ |  $5.4$%
lmart |  $7.19$ $\pm$ |  $4.2$% |  $14$ / | $100$ |  $167$ $\pm$ |  $9.0$% |  $59.98$ $\pm$ |  $3.9$% |  $0$ / | $100$ |  $1925$ $\pm$ |  $4.9$% |  $2165.96$ $\pm$ |  $2.0$% |  $0$ / | $54$ |  $45\,319$ $\pm$ |  $3.4$%
gcnn |  $6.59$ $\pm$ |  $3.1$% |  $85$ / | $100$ |  $134$ $\pm$ |  $7.6$% |  $42.48$ $\pm$ |  $2.7$% |  $100$ / | $100$ |  $1450$ $\pm$ |  $3.3$% |  $1489.91$ $\pm$ |  $3.3$% |  $66$ / | $70$ |  $29\,981$ $\pm$ |  $4.9$%
| Set Covering
fsb |  $4.11$ $\pm$ |  $12.1$% |  $0$ / | $100$ |  $6$ $\pm$ |  $30.3$% |  $86.90$ $\pm$ |  $12.9$% |  $0$ / | $100$ |  $72$ $\pm$ |  $19.4$% |  $1813.33$ $\pm$ |  $5.1$% |  $0$ / | $68$ |  $400$ $\pm$ |  $7.5$%
rpb |  $2.74$ $\pm$ |  $7.8$% |  $0$ / | $100$ |  $10$ $\pm$ |  $32.1$% |  $17.41$ $\pm$ |  $6.6$% |  $0$ / | $100$ |  $689$ $\pm$ |  $21.2$% |  $136.17$ $\pm$ |  $7.9$% |  $13$ / | $100$ |  $5511$ $\pm$ |  $11.7$%
trees |  $2.47$ $\pm$ |  $7.3$% |  $0$ / | $100$ |  $86$ $\pm$ |  $15.9$% |  $23.70$ $\pm$ |  $11.2$% |  $0$ / | $100$ |  $976$ $\pm$ |  $14.4$% |  $451.39$ $\pm$ |  $14.6$% |  $0$ / | $95$ |  $10\,290$ $\pm$ |  $16.2$%
svmrank |  $2.31$ $\pm$ |  $6.8$% |  $0$ / | $100$ |  $77$ $\pm$ |  $15.0$% |  $23.10$ $\pm$ |  $9.8$% |  $0$ / | $100$ |  $867$ $\pm$ |  $13.4$% |  $364.48$ $\pm$ |  $7.7$% |  $0$ / | $98$ |  $6329$ $\pm$ |  $7.7$%
lmart |  $1.79$ $\pm$ |  $6.0$% |  $75$ / | $100$ |  $77$ $\pm$ |  $14.9$% |  $14.42$ $\pm$ |  $9.5$% |  $1$ / | $100$ |  $873$ $\pm$ |  $14.3$% |  $222.54$ $\pm$ |  $8.6$% |  $0$ / | $100$ |  $7006$ $\pm$ |  $6.9$%
gcnn |  $1.85$ $\pm$ |  $5.0$% |  $25$ / | $100$ |  $70$ $\pm$ |  $12.0$% |  $10.29$ $\pm$ |  $7.1$% |  $99$ / | $100$ |  $657$ $\pm$ |  $12.2$% |  $114.16$ $\pm$ |  $10.3$% |  $87$ / | $100$ |  $5169$ $\pm$ |  $14.9$%
| Combinatorial Auction
fsb |  $30.36$ $\pm$ |  $19.6$% |  $4$ / | $100$ |  $14$ $\pm$ |  $34.5$% |  $214.25$ $\pm$ |  $15.2$% |  $1$ / | $100$ |  $76$ $\pm$ |  $15.8$% |  $742.91$ $\pm$ |  $9.1$% |  $15$ / | $90$ |  $55$ $\pm$ |  $7.2$%
rpb |  $26.55$ $\pm$ |  $16.2$% |  $9$ / | $100$ |  $22$ $\pm$ |  $31.9$% |  $156.12$ $\pm$ |  $11.5$% |  $8$ / | $100$ |  $142$ $\pm$ |  $20.6$% |  $631.50$ $\pm$ |  $8.1$% |  $14$ / | $96$ |  $110$ $\pm$ |  $15.5$%
trees |  $28.96$ $\pm$ |  $14.7$% |  $3$ / | $100$ |  $135$ $\pm$ |  $20.0$% |  $159.86$ $\pm$ |  $15.3$% |  $3$ / | $100$ |  $401$ $\pm$ |  $11.6$% |  $671.01$ $\pm$ |  $11.1$% |  $1$ / | $95$ |  $381$ $\pm$ |  $11.1$%
svmrank |  $23.58$ $\pm$ |  $14.1$% |  $11$ / | $100$ |  $117$ $\pm$ |  $20.5$% |  $130.86$ $\pm$ |  $13.6$% |  $13$ / | $100$ |  $348$ $\pm$ |  $11.4$% |  $586.13$ $\pm$ |  $10.0$% |  $21$ / | $95$ |  $321$ $\pm$ |  $8.8$%
lmart |  $23.34$ $\pm$ |  $13.6$% |  $16$ / | $100$ |  $117$ $\pm$ |  $20.7$% |  $128.48$ $\pm$ |  $15.4$% |  $23$ / | $100$ |  $349$ $\pm$ |  $12.9$% |  $582.38$ $\pm$ |  $10.5$% |  $15$ / | $95$ |  $314$ $\pm$ |  $7.0$%
gcnn |  $22.10$ $\pm$ |  $15.8$% |  $57$ / | $100$ |  $107$ $\pm$ |  $21.4$% |  $120.94$ $\pm$ |  $14.2$% |  $52$ / | $100$ |  $339$ $\pm$ |  $11.8$% |  $563.36$ $\pm$ |  $10.7$% |  $30$ / | $95$ |  $338$ $\pm$ |  $10.9$%
| Capacitated Facility Location
fsb |  $23.58$ $\pm$ |  $29.9$% |  $9$ / | $100$ |  $7$ $\pm$ |  $35.9$% |  $1503.55$ $\pm$ |  $20.9$% |  $0$ / | $74$ |  $38$ $\pm$ |  $28.2$% |  $3600.00$ $\pm$ |  $0.0$% |  $0$ / | $0$ |  n/a $\pm$ | n/a%
rpb |  $8.77$ $\pm$ |  $11.8$% |  $7$ / | $100$ |  $20$ $\pm$ |  $36.1$% |  $110.99$ $\pm$ |  $24.4$% |  $41$ / | $100$ |  $729$ $\pm$ |  $37.3$% |  $2045.61$ $\pm$ |  $18.3$% |  $22$ / | $42$ |  $2675$ $\pm$ |  $24.0$%
trees |  $10.75$ $\pm$ |  $22.1$% |  $1$ / | $100$ |  $76$ $\pm$ |  $44.2$% |  $1183.37$ $\pm$ |  $34.2$% |  $1$ / | $47$ |  $4664$ $\pm$ |  $45.8$% |  $3565.12$ $\pm$ |  $1.2$% |  $0$ / | $3$ |  $38\,296$ $\pm$ |  $4.1$%
svmrank |  $8.83$ $\pm$ |  $14.9$% |  $2$ / | $100$ |  $46$ $\pm$ |  $32.2$% |  $242.91$ $\pm$ |  $29.3$% |  $1$ / | $96$ |  $546$ $\pm$ |  $26.0$% |  $2902.94$ $\pm$ |  $9.6$% |  $1$ / | $18$ |  $6256$ $\pm$ |  $15.1$%
lmart |  $7.31$ $\pm$ |  $12.7$% |  $30$ / | $100$ |  $52$ $\pm$ |  $38.1$% |  $219.22$ $\pm$ |  $36.0$% |  $15$ / | $91$ |  $747$ $\pm$ |  $35.1$% |  $3044.94$ $\pm$ |  $7.0$% |  $0$ / | $12$ |  $8893$ $\pm$ |  $3.5$%
gcnn |  $6.43$ $\pm$ |  $11.6$% |  $51$ / | $100$ |  $43$ $\pm$ |  $40.2$% |  $192.91$ $\pm$ |  $110.2$% |  $42$ / | $82$ |  $1841$ $\pm$ |  $88.0$% |  $2024.37$ $\pm$ |  $30.6$% |  $25$ / | $29$ |  $2997$ $\pm$ |  $26.3$%
| Maximum Independent Set

#### Training

We train the models on each benchmark separately. Namely, for each benchmark,
we generate 100,000 branching samples extracted from 10,000 randomly generated
instances for training, 20,000 branching samples from 2,000 instances for
validation, and same for test (see supplementary materials for details). We
report in Table 1 the test accuracy of each machine learning model over five
seeds, as the percentage of times the highest ranked decision of the model
(acc@1), one of the five highest (acc@5) and one of the ten highest (acc@10)
is a variable which is given the highest strong branching score.

#### Evaluation

Evaluation is performed for each problem difficulty (Easy, Medium, Hard) on 20
new instances using five different seeds111In addition to ML models which are
re-trained with a different seed, all major MILP solvers have a parameter,
_seed_ , that randomizes some tie-breaking rules, so as to be able to report
aggregated results over the same instance., which amounts to a total of 100
solving attempts per method. We report standard metrics for MILP
benchmarking222See e.g. http://plato.asu.edu/bench.html, that is: the
1-shifted geometric mean of the solving times in seconds, including running
times of unsolved instances without extra penalization (Time); the hardware-
independent final node counts of instances that are solved by all baselines
(Nodes); and the number of times each branching policy results in the fastest
solving time, over the number of instances solved (Win). Policy evaluation
results are displayed in Table 2. Note that we also report the average per-
instance standard deviation, so “$64\pm 13.6\%$ nodes” means it took on
average 64 nodes to solve an instance, and when solving one of those instances
the number of nodes varied by 13.6% on average.

###  5.2 Comparative experiment

In terms of prediction accuracy (Table 1), gcnn clearly outperforms the
baseline competitors on all four problems, while svmrank and lmart are on par
with each other and the performance of trees is the lowest.

At solving time (Table 2), the accuracy of each method is clearly reflected in
the number of nodes required to solve the instances. Interestingly, the best
method in terms of nodes is not necessarily the best in terms of total solving
time, which also takes into account the computational cost of each branching
policy, i.e., the feature extraction and inference time. The svmrank approach,
despite being slightly better than lmart in terms of number of nodes, is also
slower due to a worse running time / number of nodes trade-off. Our gcnn model
clearly dominates overall, except on combinatorial auction (Easy) and maximum
independent set (Medium) instances, where lmart and rpb are respectively
faster.

Our gcnn model generalizes well to instances of size larger than seen during
training, and outperforms SCIP’s default branching rule rpb in terms of
running time in almost every configuration. In particular and strikingly, it
significantly outperforms rpb in terms of nodes on medium and hard instances
for setcover and combinatorial auction problems. As expected, the fsb expert
brancher is not competitive in terms of running time, despite producing very
small search trees. The maximum independent set problem seems particularly
challenging for generalization, as all machine learning approaches report a
lower number of solved instances than the default rpb brancher, and gcnn,
despite being the fastest machine learning approach overall, exhibits a high
variability both in terms of time and number of nodes.

For the first time in the literature a machine-learning-based approach is
compared with an essentially full-fledged MILP solver. For this reason, the
results are particularly impressive, and indicate that gcnn is a very serious
candidate to be implemented within a MILP solver, as an additional tool to
speed up mixed-integer linear programming solvers. Also, they suggest that
more could be gained from a tight integration within a complex software, like
any MILP solver is.

###  5.3 Ablation study

We present an ablation study of our proposed GCNN model on the set covering
problem by comparing three variants of the convolution operation in (4): mean
rather than sum convolutions (mean), sum convolutions without our prenorm
layer (sum) and finally sum convolutions with prenorm layers, which is the
model we use throughout our experiments (gcnn).

Results on test instances are reported in Table 3. The solving performance of
both variants mean and sum is very similar to that of our baseline gcnn on
small instances. On large instances however, the variants perform
significantly worse in terms of both solving time and number of nodes,
especially on hard instances. This empirical evidence supports our hypothesis
that sum-convolutions offer a better architectural prior than mean-convolution
from the task of learning to branch, and that our prenorm layer helps for
stabilizing training and improving generalization.

Table 3: Ablation study of our GCNN model on the set covering problem. Sum convolutions generalize better to larger instances, especially when combined with a prenorm layer. | Accuracies | Easy | Medium | Hard
---|---|---|---|---
Model | acc@1 | acc@5 | acc@10 | time | wins | nodes | time | wins | nodes | time | wins | nodes
mean |  $65.4$ $\pm$ | $0.1$ |  $92.4$ $\pm$ | $0.1$ |  $98.2$ $\pm$ | $0.0$ |  $6.7$ $\pm$ |  $3$% |  $13$ / | $100$ |  $134$ $\pm$ |  $6$% |  $43.7$ $\pm$ |  $3$% |  $19$ / | $100$ |  $1894$ $\pm$ |  $4$% |  $1593.0$ $\pm$ |  $4$% |  $6$ / | $70$ |  $62\,227$ $\pm$ |  $6$%
sum |  $65.5$ $\pm$ | $0.2$ |  $92.3$ $\pm$ | $0.2$ |  $98.1$ $\pm$ | $0.1$ |  $6.6$ $\pm$ |  $3$% |  $27$ / | $100$ |  $134$ $\pm$ |  $6$% |  $42.5$ $\pm$ |  $3$% |  $45$ / | $100$ |  $1882$ $\pm$ |  $4$% |  $1511.7$ $\pm$ |  $3$% |  $22$ / | $70$ |  $57\,864$ $\pm$ |  $4$%
gcnn |  $65.5$ $\pm$ | $0.1$ |  $92.4$ $\pm$ | $0.1$ |  $98.2$ $\pm$ | $0.0$ |  $6.6$ $\pm$ |  $3$% |  $60$ / | $100$ |  $134$ $\pm$ |  $8$% |  $42.5$ $\pm$ |  $3$% |  $36$ / | $100$ |  $1870$ $\pm$ |  $3$% |  $1489.9$ $\pm$ |  $3$% |  $42$ / | $70$ |  $56\,348$ $\pm$ |  $5$%

##  6 Discussion

The objective of branch-and-bound is to solve combinatorial optimization
problems as fast as possible. Branching policies must therefore balance the
quality of decisions taken with the time spent to take each decision. An
extreme example of this tradeoff is strong branching: this policy takes
excellent decisions leading to low number of nodes overall, but every
decision-making step is so slow that the overall running time is not
competitive. Early experiments showed that we could take better decisions and
decrease the number of nodes slightly on average by training a GCNN policy
with more layers or with a larger embedding size. However, this would also
lead to increased computational costs for inference and slightly larger times
at each decision, and in the end increased solving times on average. The
policy architecture we propose is thus a compromise between learning capacity
and inference speed, something that is not traditionally a concern within the
machine learning community.

Another concern among the combinatorial optimization community is the ability
of policies trained on small instances to generalize to larger instances. We
were able to show that machine learning methods, and the GCNN model in
particular, can generalize to fairly larger instances. However, in general it
is expected that the improvement in performance decreases as our model is
evaluated on progressively larger problems, as can already be observed from
Table 2. In early experiments with even larger instances (huge), we observed a
performance drop for the model trained on our small instances. This could
presumably be remedied by training on larger instances in the first place, and
indeed a model trained on medium instances did perform well on those huge
instances again. In any case, there are limits as to the generalization
ability of any learned branching policy, and since the limit is likely very
dependent on the problem structure, it is difficult to give any precise
quantitative estimates a priori. This desirable ability to generalize outside
of the training distribution, sometimes termed transfer learning, is also not
a traditional concern in the machine learning community.

##  7 Conclusion

We formulated branch-and-bound, the standard exact method for solving mixed-
integer linear programs, as a Markov decision process. In this context, we
proposed and evaluated a novel approach for tackling the branching problem, by
expressing the state of the branch-and-bound process as a bipartite graph,
which reduces the need for feature engineering by naturally leveraging the
variable-constraint structure of MILP problems, and allows for the encoding of
branching policies as a graph convolutional neural network. We demonstrated on
four NP-hard problems that, by adopting a simple imitation learning scheme,
the policies learned by a GCNN were outperforming previously proposed machine
learning approaches for branching, and could also outperform the default
branching strategy of SCIP, a modern open-source solver. Most importantly, we
demonstrated that the learned policies could generalize to instance sizes
larger than seen during training. This is essential since collecting strong
branching decisions, hence training, can be computationally prohibitive on
large instances. Our work indicates that the GCNN model, especially using sum
convolutions with the proposed prenorm layer, is a good architectural prior
for the task of branching in MILP.

In future work, we would like to assess the viability of our approach on a
broader set on combinatorial problems, and also experiment with reinforcement
learning methods for improving over the policies learned by imitation. Also,
we believe that there is plenty of room for hydrid approaches combining
traditional methods and machine learning for branching, and we would like to
dig deeper into the learned policies in order to extract some knowledge of
interest for the MILP community.

## Acknowledgements

We would like to thank the anonymous reviewers whose contributions helped
considerably improve the quality of this paper. We would also like to thank
Ambros Gleixner and Benjamin Müller for enlightening discussions and technical
help regarding SCIP, as well as Gonzalo Muñoz, Aleksandr Kazachkov and Giulia
Zarpellon for insightful discussions on variable selection. Finally, we thank
Jason Jo, Meng Qu and Mike Pieper for their helpful comments on the structure
of the paper.

This work was supported by the Canada First Research Excellence Fund (CFREF),
IVADO, CIFAR, GERAD, and Canada Excellence Research Chairs (CERC).

## References

  * Achterberg and Berthold [2009] Tobias Achterberg and Timo Berthold.  Hybrid branching.  In _Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems_ , 2009.
  * Achterberg and Wunderling [2013] Tobias Achterberg and Roland Wunderling.  Mixed integer programming: Analyzing 12 years of progress.  _Facets of Combinatorial Optimization_ , pages 449–481, 2013.
  * Achterberg et al. [2005] Tobias Achterberg, Thorsten Koch, and Alexander Martin.  Branching rules revisited.  _Operations Research Letters_ , 33:42–54, 2005.
  * Alvarez et al. [2017] Alejandro M. Alvarez, Quentin Louveaux, and Louis Wehenkel.  A machine learning-based approximation of strong branching.  _INFORMS Journal on Computing_ , 29:185–195, 2017.
  * Applegate et al. [1995] David Applegate, Robert Bixby, Vašek Chvátal, and William Cook.  Finding cuts in the TSP.  Technical report, DIMACS, 1995.
  * Åström [1965] Karl J. Åström.  Optimal control of Markov processes with incomplete state information.  _Journal of Mathematical Analysis and Applications_ , 10:174–205, 1965.
  * Balas and Ho [1980] Egon Balas and Andrew Ho.  Set covering algorithms using cutting planes, heuristics, and subgradient optimization: a computational study.  In _Combinatorial Optimization_ , pages 37–60. Springer, 1980.
  * Balcan et al. [2018] Maria-Florina Balcan, Travis Dick, Tuomas Sandholm, and Ellen Vitercik.  Learning to branch.  In _Proceedings of the International Conference on Machine Learning_ , 2018.
  * Balunovic et al. [2018] Mislav Balunovic, Pavol Bielik, and Martin Vechev.  Learning to solve SMT formulas.  In _Advances in Neural Information Processing Systems 31_ , pages 10338–10349, 2018.
  * Bengio et al. [2018] Yoshua Bengio, Andrea Lodi, and Antoine Prouvost.  Machine learning for combinatorial optimization: a methodological tour d’horizon.  _arXiv:1811.06128_ , 2018.
  * Bergman et al. [2016] David Bergman, Andre A. Cire, Willem-Jan Van Hoeve, and John Hooker.  _Decision diagrams for optimization_.  Springer, 2016.
  * Bruna et al. [2014] Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun.  Spectral networks and locally connected networks on graphs.  In _Proceedings of the Second International Conference on Learning Representations_ , 2014.
  * Burges [2010] Christopher J. C. Burges.  From RankNet to LambdaRank to LambdaMART: An Overview.  Technical report, Microsoft Research, 2010.
  * Cornuejols et al. [1991] Gerard Cornuejols, Ramaswami Sridharan, and Jean Michel Thizy.  A comparison of heuristics and relaxations for the capacitated plant location problem.  _European Journal of Operational Research_ , 50:280 – 297, 1991.
  * Di Liberto et al. [2016] Giovanni Di Liberto, Serdar Kadioglu, Kevin Leo, and Yuri Malitsky.  Dash: Dynamic approach for switching heuristics.  _European Journal of Operational Research_ , 248:943–953, 2016.
  * Duvenaud et al. [2015] David K. Duvenaud, Dougal Maclaurin, Jorge Aguilera-Iparraguirre, Rafael Bombarell, Timothy Hirzel, Alán Aspuru-Guzik, and Ryan P. Adams.  Convolutional networks on graphs for learning molecular fingerprints.  In _Advances in Neural Information Processing Systems 28_ , pages 2224–2232, 2015.
  * Fischetti and Monaci [2012] Matteo Fischetti and Michele Monaci.  Branching on nonchimerical fractionalities.  _Operations Research Letters_ , 40:159 – 164, 2012.
  * Geurts et al. [2006] Pierre Geurts, Damien Ernst, and Louis Wehenkel.  Extremely randomized trees.  _Machine learning_ , 63:3–42, 2006.
  * Gilmer et al. [2017] Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl.  Neural message passing for quantum chemistry.  In _Proceedings of the Thirty-Fourth International Conference on Machine Learning_ , pages 1263–1272, 2017.
  * Gleixner et al. [2018] Ambros Gleixner, Michael Bastubbe, Leon Eifler, Tristan Gally, Gerald Gamrath, Robert Lion Gottwald, Gregor Hendel, Christopher Hojny, Thorsten Koch, Marco E. Lübbecke, Stephen J. Maher, Matthias Miltenberger, Benjamin Müller, Marc E. Pfetsch, Christian Puchert, Daniel Rehfeldt, Franziska Schlösser, Christoph Schubert, Felipe Serrano, Yuji Shinano, Jan Merlin Viernickel, Matthias Walter, Fabian Wegscheider, Jonas T. Witt, and Jakob Witzig.  The SCIP Optimization Suite 6.0.  Zib-report, Zuse Institute Berlin, July 2018.
  * Gleixner et al. [2019] Ambros Gleixner, Gregor Hendel, Gerald Gamrath, Tobias Achterberg, Michael Bastubbe, Timo Berthold, Philipp Christophel, Kati Jarck, Thorsten Koch, Jeff Linderoth, Marco Lübbecke, Hans D. Mittelmann, Derya Ozyurt, Ted K. Ralphs, Domenico Salvagnin, and Yuji Shinano.  MIPLIB 2017: Data-Driven Compilation of the 6th Mixed-Integer Programming Library.  Technical report, Optimization Online, August 2019.
  * Glorot and Bengio [2010] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , pages 249–256, 2010.
  * Gori et al. [2005] Marco Gori, Gabriele Monfardini, and Franco Scarselli.  A new model for learning in graph domains.  In _Proceedings of the 2005 IEEE International Joint Conference on Neural Networks_ , volume 2, pages 729–734, 2005.
  * Hansknecht et al. [2018] Christoph Hansknecht, Imke Joormann, and Sebastian Stiller.  Cuts, primal heuristics, and learning to branch for the time-dependent traveling salesman problem.  arXiv:1805.01415, 2018.
  * He et al. [2014] He He, Hal III Daumé, and Jason Eisner.  Learning to search in branch-and-bound algorithms.  In _Advances in Neural Information Processing Systems 27_ , pages 3293–3301, 2014.
  * Howard [1960] Ronald A. Howard.  _Dynamic Programming and Markov Processes_.  MIT Press, Cambridge, MA, 1960.
  * Hussein et al. [2017] Ahmed Hussein, Mohamed M. Gaber, Eyad Elyan, and Chrisina Jayne.  Imitation learning: A survey of learning methods.  _ACM Computing Surveys_ , 50:21, 2017.
  * Joachims [2002] Thorsten Joachims.  Optimizing search engines using clickthrough data.  In _Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 133–142, 2002.
  * Karzan et al. [2009] Fatma K. Karzan, George L. Nemhauser, and Martin W. P. Savelsbergh.  Information-based branching schemes for binary linear mixed integer problems.  _Mathematical Programming Computation_ , 1:249–293, 2009\.
  * Khalil et al. [2016] Elias B. Khalil, Pierre Le Bodic, Le Song, George Nemhauser, and Bistra Dilkina.  Learning to branch in mixed integer programming.  In _Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence_ , pages 724–731, 2016.
  * Khalil et al. [2017a] Elias B. Khalil, Hanjun Dai, Yuyu Zhang, Bistra Dilkina, and Le Song.  Learning combinatorial optimization algorithms over graphs.  In _Advances in Neural Information Processing Systems 30_ , pages 6348–6358, 2017a.
  * Khalil et al. [2017b] Elias B. Khalil, Bistra Dilkina, George L. Nemhauser, Shabbir Ahmed, and Yufen Shao.  Learning to run heuristics in tree search.  In _Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence_ , pages 659–666, 2017b.
  * Kipf and Welling [2017] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In _Proceedings of the Fifth International Conference on Learning Representations_ , 2017.
  * Kulesza and Pereira [2008] Alex Kulesza and Fernando Pereira.  Structured learning with approximate inference.  In _Advances in Neural Information Processing Systems 21_ , pages 785–792, 2008.
  * Land and Doig [1960] Ailsa H. Land and Alison G. Doig.  An automatic method of solving discrete programming problems.  _Econometrica_ , 28:497–520, 1960.
  * Lederman et al. [2018] Gil Lederman, Markus N. Rabe, and Sanjit A. Seshia.  Learning heuristics for automated reasoning through deep reinforcement learning.  arXiv:1807.08058, 2018.
  * Leyton-Brown et al. [2000] Kevin Leyton-Brown, Mark Pearson, and Yoav Shoham.  Towards a universal test suite for combinatorial auction algorithms.  In _Proceedings of the Second ACM Conference on Electronic Commerce_ , pages 66–76, 2000.
  * Li et al. [2016] Yujia Li, Daniel Tarlow, Marc Brockschmidt, and Richard S. Zemel.  Gated graph sequence neural networks.  In Yoshua Bengio and Yann LeCun, editors, _Proceedings of the Fourth International Conference on Learning Representations_ , 2016.
  * Liang et al. [2016] Jia Hui Liang, Vijay Ganesh, Pascal Poupart, and Krzysztof Czarnecki.  Learning rate based branching heuristic for SAT solvers.  In _International Conference on Theory and Applications of Satisfiability Testing_ , pages 123–140, 2016.
  * Linderoth and Savelsbergh [1999] Jeff T. Linderoth and Martin W. P. Savelsbergh.  A computational study of search strategies for mixed integer programming.  _INFORMS Journal on Computing_ , 11:173–187, 1999.
  * Lodi and Zarpellon [2017] Andrea Lodi and Giulia Zarpellon.  On learning and branching: a survey.  _TOP_ , 25:207–236, 2017.
  * Paschos [2014] Vangelis T. Paschos, editor.  _Applications of Combinatorial Optimization_.  Mathematics and Statistics. Wiley-ISTE, second edition, 2014.
  * Patel and Chinneck [2007] Jagat Patel and John W. Chinneck.  Active-constraint variable ordering for faster feasibility of mixed integer linear programs.  _Mathematical Programming_ , 110:445–474, 2007.
  * Pochet and Wolsey [2006] Yves Pochet and Laurence A. Wolsey.  _Production planning by mixed integer programming_.  Springer Science and Business Media, 2006.
  * Pomerleau [1991] Dean A. Pomerleau.  Efficient training of artificial neural networks for autonomous navigation.  _Neural Computation_ , 3:88–97, 1991.
  * Scarselli et al. [2009] Franco Scarselli, Marco Gori, Ah C. Tsoi, Markus Hagenbuchner, and Gabriele Monfardini.  The graph neural network model.  _IEEE Trans. Neural Networks_ , 20(1):61–80, 2009\.
  * Selsam et al. [2019] Daniel Selsam, Matthew Lamm, Benedikt Bünz, Percy Liang, Leonardo de Moura, and David L. Dill.  Learning a sat solver from single-bit supervision.  In _Proceedings of the Seventh International Conference on Learning Representations_ , 2019.
  * Song et al. [2018] Jialin Song, Ravi Lanka, Albert Zhao, Yisong Yue, and Masahiro Ono.  Learning to search via retrospective imitation.  arXiv:1804.00846, 2018.
  * Srikumar and Roth [2012] Vivek G. K. Srikumar and Dan Roth.  On amortizing inference cost for structured prediction.  In _Proceedings of the 2012 Joint Conference on Empirical Methods in Natural Language Processing and Computational Natural Language Learning_ , pages 1114–1124, 2012.
  * Sutton and Barto [2018] Richard S. Sutton and Andrew G. Barto.  _Reinforcement Learning: An Introduction_.  MIT Press, second edition, 2018.
  * Wainwright et al. [2005] Martin J. Wainwright, Tommi S. Jaakkola, and Alan S. Willsky.  MAP estimation via agreement on trees: message-passing and linear programming.  _IEEE Transactions on Information Theory_ , 51:3697–3717, 2005.
  * Wolsey [1988] Laurence A. Wolsey.  _Integer Programming_.  Wiley-Blackwell, 1988.
  * Xu et al. [2019] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  In _Proceedings of the Seventh International Conference on Learning Representations_ , 2019.


# Insight from the elliptic flow of identified hadrons measured in
relativistic heavy-ion collisions

Praphupada Dixit and Md. Nasim  Department of Physics, Indian Institute of
Science Education and Research, Berhampur 760010, India

###### Abstract

In this paper, we have discussed the transverse momentum dependence of
elliptic flow measured in relativistic heavy-ion collisions. Under the
assumption of quark recombination model, $v_{2}$ of light and strange quarks
are obtained from the measured $v_{2}$ of identified hadrons. The $v_{2}$ of
strange quarks obtained from $\phi$-meson $v_{2}$ is found to be similar to
that of from $\Omega$ $v_{2}$ in Au+Au collisions at $\sqrt{s_{NN}}$ = 200
GeV. This may indicates that both $\phi$ and $\Omega$ are produced through
quark recombination at top RHIC energy. However, we find the $v_{2}$ of light
quarks is consistently higher than that of strange quarks at the low
transverse momentum. We have shown, $v_{2}$ of light and strange quarks are
consistent with each other when plotted as a function of its transverse
kinetic energy. The observed scaling of $v_{2}$, works at both RHIC and LHC
energies. It is also shown that the observed scaling holds for other flow
harmonics, e.g. $v_{3}$. The new observed scaling is purely empirical and
needs more theoretical input to understand the physics behind it.

###### pacs:

25.75.Ld

##  I Introduction

Quarks are the fundamental particles or building blocks of our universe. The
bound states of quarks are known as hadrons. Due to strong coupling through
gluon exchange, quarks are always found to be bound within hadrons. No free
quarks exist in nature. The theory of quantum chromodynamics (QCD) predicts
that at very high temperature and/or high density, hadrons might get melted
into a deconfined matter of quarks and gluons. The deconfined medium of quarks
and gluons is known as quark-gluon plasma (QGP). The relativistic heavy-ion
collisions provide a unique opportunity to study the QGP in laboratory
experiments white . The medium created in the heavy-ion collision is very hot
($\sim$4 trillion degrees Celsius) and dense and also extremely short-lived
($\sim$10-23 sec). In experiments, we are only able to detect the freely
streaming final state particles emerging from the collisions. Using the
information carried by these particles as probes, we can study the properties
of the medium created in the collision. The properties of the hot and dense
medium formed in the heavy-ion collisions can be studied by using azimuthal
anisotropies of produced hadrons in momentum space. These momentum
anisotropies arise due to initial spatial anisotropies followed by multiple
partonic interactions. The azimuthal distribution of the produced particles
can be described by using a Fourier series. Then, the Fourier coefficients of
different harmonics quantify the degree of final state anisotropies. Within a
hydrodynamical framework, the $v_{2}$ has been shown to be sensitive to the
equation of the state of the system formed in the collisions flow1 ; flow2 ;
hydro_flow .

The elliptic flow parameter is defined as the $2^{\mathrm{nd}}$ Fourier
coefficient, $v_{2}$, of the particle distributions in emission azimuthal
angle ($\phi$) with respect to the reaction plane angle ($\Psi$) art , and can
be written as

$\frac{dN}{d(\phi-\Psi)}\propto 1+2v_{2}\cos(2(\phi-\Psi)).$ (1)

For a given rapidity window, the second coefficient is

$v_{2}=\langle\cos(2(\phi-\Psi))\rangle.$ (2)

The angular bracket represents the average over all the produced particles.
Experimentally $v_{2}$ of different identified hadrons are measured as a
function of transverse momentum $p_{T}$, centrality, and rapidity pidv2_bes1 ;
pidprl_bes1 ; pidprl_200 ; rhicflow ; phobosflow ; starflow ; phenixflow .

In this paper, we have discussed $v_{2}(p_{T})$ of identified hadrons measured
in heavy-ion collisions and its physics implications.

##  II Results

Experimentally, $v_{2}$ of different identified hadrons are measured as a
function of transverse momentum at Relativistic Heavy-Ion Collider (RHIC) and
Large Hadron Collider (LHC). The measurements show a mass ordering in $v_{2}$
for $p_{T}$ $<$ 2 GeV/c and particle type dependence (i.e baryon vs. meson) at
the intermediate $p_{T}$ (2.0 $<$ $p_{T}$ $<$ 4.0 GeV/c). The observed mass
ordering is the effect of radial flow of the medium created in the collision
flow2 ; hydro_flow  and particle-type dependence can be explained by the
quarks coalescence/recombination models coal_fries ; coal_greco ;
coal_voloshin ; coal_all_1 ; coal_all_2 ; coal_all_3 ; coal_all_4 .

Within the framework of the coalescence mechanism, the $v_{2}$ of hadrons and
quarks can be related by the following equation:

$v_{2}(p_{T})=n_{q}v_{2}^{q}(p_{T}/n_{q}),$ (3)

where $p_{T}$ is the transverse momentum of hadrons, $n_{q}$ is the number of
constituent quarks coal_fries . The $v_{2}$ and $v_{2}^{q}$ are the elliptic
flow of hadrons and quarks, respectively.

Figure 1: (Color online) The $v_{2}$ of light and strange quarks obtained from
proton, $\phi$, $\Omega$, $\Lambda$, and $\Xi$ as a function of $p_{T}/n_{q}$
in Au+Au collisions at $\sqrt{s_{NN}}$ = 200 GeV pidprl_200 .

Figure 1 shows the $n_{q}$ scaled $v_{2}$ as a function of $p_{T}/n_{q}$ in
Au+Au collisions at $\sqrt{s_{NN}}$ = 200 GeV coal_review . We observe that
$n_{q}$ scaled $\phi$ $v_{2}$ is similar to $n_{q}$ scaled $\Omega$ $v_{2}$
within the measurement uncertainty. According to Ref. phi_omega_star  the
production of $\phi$ and $\Omega$ is mainly through the coalescence of strange
quarks at these energies, therefore both $n_{q}$ scaled $\phi$ and $\Omega$
$v_{2}$ should give $v_{2}$ of strange quarks and they should match. Our
observation is consistent with the picture of $\phi$ and $\Omega$ production
through quark coalescence.

Now, $n_{q}$ scaled proton $v_{2}$ should be equivalent to $v_{2}$ of light
quarks (considering up and down quarks have the same mass). We have also
calculated $v_{2}$ of light quarks from the measured $v_{2}$ of $\Xi$ and
$\Lambda$ using the following equation coal_all_4 ; coal_all_5 ,

$v_{2}^{u}(p_{T})=\frac{1}{3}[2v_{2}^{\Lambda}((2+r)p_{T})-v_{2}^{\Xi}((1+2r)p_{T})];$
(4)

where the factor $r=m_{s}/m_{u}$. Using constituent quark masses, the factor
$r$ becomes 1.667. We can see from Fig. 1, the light quarks $v_{2}$ at 200 GeV
obtained using Eq. 4 is equivalent to that obtained from the proton $v_{2}$.
From Fig. 1, we see $n_{q}$ scaled proton $v_{2}$ is higher than $n_{q}$
scaled $\phi$ and $\Omega$ $v_{2}$ at low $p_{T}$. This could mean $v_{2}$ of
strange quarks is smaller than $v_{2}$ of light quarks at low $p_{T}$ or this
could be due to small hadronic interaction cross-sections of $\phi$ and
$\Omega$ compared to non-strange hadrons pidprl_200 ; refree_1 ; refree_3 .
Alternatively, the difference between $n_{q}$ scaled proton $v_{2}$ and
$n_{q}$ scaled $\phi$ $v_{2}$ at low $p_{T}$ could be due to the effect of
radial flow on quarks $v_{2}$ at the partonic stage.

Figure 2: (Color online) Panel (a) : The $v_{2}$ of light and strange quarks
obtained from proton, $\phi$, $\Omega$, $\Lambda$, and $\Xi$ as a function of
$p_{T}/n_{q}$ in Au+Au collisions at $\sqrt{s_{NN}}$ = 200 GeV. Panel (b) :
The $v_{2}$ of quarks as a function of $KE_{T}^{q}$=
$\sqrt{(p_{T}/n_{q})^{2}+m_{q}^{2}}-m_{q}$, where $p_{T}$ is the transverse
momentum of hadron. Here we consider $m_{q}$ = 4 MeV/c2 and 100 MeV/c2 for
light and strange quarks, respectively. Panel (c) : Same as panel (b), but
here $m_{q}$ = 140 MeV/c2 for strange quarks. Blue dashed curves are the fit
to light quarks $v_{2}$ using Eq. 5 and red dashed lines are the constant
polynomial fit to the ratios shown in the respective bottom panels. Figure 3:
(Color online) $\chi^{2}/ndf$ for different s-quark mass for the best match
between strange and light quarks $v_{2}$ as a function of transverse kinetic
energy. This is for Au+Au collisions at $\sqrt{s_{NN}}$ =200 GeV. Panel (a):
s-quarks $v_{2}$ derived from the $\Omega$ is used. Panel (b): s-quarks
$v_{2}$ derived from the $\phi$ is used.

The $v_{2}$ as a function of the transverse kinetic energy of particles
produced in ultra-relativistic heavy-ion collisions is found to be a more
effective way of showing scaling among different hadrons, as shown by the
PHENIX experiment scaling_phenix . It was shown that $v_{2}/n_{q}$ as a
function of $KE_{T}$/$n_{q}$ shows better scaling compared to $v_{2}/n_{q}$ as
a function of $p_{T}/n_{q}$. Here $KE_{T}$ = $m_{T}$ -$m_{0}$ is the
transverse kinetic energy of the hadrons. However, this scaling has no strong
theory background.

In this paper, we have shown an alternative empirical scaling relation between
$v_{2}/n_{q}$ as a function of $KE_{T}^{q}$ =
$\sqrt{(p_{T}/n_{q})^{2}+m_{q}^{2}}-m_{q}$, where $p_{T}$ is the transverse
momentum of hadron. If we use $m_{q}$ as the bare mass of a quark and
$p_{T}/n_{q}$ as its transverse momentum, then the quantity $KE_{T}^{q}$ is
equivalent to transverse kinetic energy of quark. The values of bare or
current masses of up ($u$), down ($d$) and strange ($s$) quarks are 1.8-2.8
MeV/$c^{2}$, 4.3-5.2 MeV/$c^{2}$ and 92-104 MeV/$c^{2}$, respectively pdg .
Since mass of $u$ and $d$ quarks are very close, in this study we took 4
MeV/$c^{2}$ as an average value of $u$ and $d$ quarks mass. Figure. 2 (middle
panel), shows $v_{2}/n_{q}$ as a function of $KE_{T}^{q}$ using mass of
strange ($m_{s}$) quarks is equal to 100 MeV/$c^{2}$. We can see the
$v_{2}/n_{q}$ as a function of $KE_{T}^{q}$ shown better scaling compared to
$v_{2}/n_{q}$ as a function of $p_{T}/n_{q}$ (left panel). However, we have
seen the best scaling obtained when we use $m_{s}$ $\sim$ 130-140 MeV/$c^{2}$.
For different $s$-quark mass, we have calculated $\chi^{2}/ndf$ to get the
best match between strange and light quarks $v_{2}$ which is shown in Fig. 3.
Since we have s-quark $v_{2}$ derived from both $\phi$ and $\Omega$, we have
shown $\chi^{2}/ndf$ distribution for $\phi$ and $\Omega$ separately. We can
see from Fig. 3, that $s$-quark mass $m_{s}$ = 130-140 MeV/c2 gives best
$\chi^{2}/ndf$ while comparing with light quark $v_{2}$ as a function of
transverse kinetic energy. The value of s-quark mass $m_{s}$ = 130-140 MeV/c2
is purely empirical and does not yet have a theoretical background. Hence, in
this paper, we present the scaling relation using both $m_{s}$ = 100 MeV/c2
(PDG value) and 140 MeV/c2 (best $\chi^{2}/ndf$).

Further, to quantify the scaling, we have shown the ratio between light and
strange quarks $v_{2}$ in the respective bottom panels. To calculate the
ratios, we have fitted the light quarks $v_{2}$ with the following function,

$v_{2}(x)=\frac{an}{1+\exp(-(xn-b)/c)}-dn,$ (5)

where $a$, $b$, $c$, $d$ and $n$ are the fit parameters. The ratios are then
fitted with a constant polynomial to get the average deviation ($p_{0}$) of
the scaling.
In Fig. 2, the pion $v_{2}$ is not included since the $n_{q}$ scaled $v_{2}$
of pions is found to be higher than that of protons, which could be due to the
effect of resonance decay in pion as reported in Ref. coal_all_4 .

###  II.1 The test of scaling at RHIC BES energies

The validity of the observed scaling is examined at other BES energies.
Surprisingly, we find that such scaling works well even at different energies.
Figure 4 shows $v_{2}/n_{q}$ as a function of $KE_{T}^{q}$ in 0-80% minimum
bias Au+Au collisions at $\sqrt{s_{NN}}$ = 39, 19.6, and 11.5 GeV. At lower
beam energies, we have used anti-proton $v_{2}$ instead of proton to avoid the
effect of transported quarks. We have seen the new proposed scaling works
equally good for both 39 and 19.6 GeV, for 11.5 GeV we need high statistics
measurements to make any conclusion. The upcoming new data from RHIC Beam
Energy Scan Phase-II (BES-II) will be useful to check this empirical scaling
at lower energies, $\sqrt{s_{NN}}$ $<$19.6 GeV.

Figure 4: (Color online) The left panel : The $v_{2}$ of light and strange
quarks obtained from $\bar{p}$ and $\phi$ (and $\Omega$, $\bar{\Lambda}$, and
$\bar{\Xi}$ for 39 GeV only) as a function of $p_{T}/n_{q}$ in 0-80% central
Au+Au collisions at $\sqrt{s_{NN}}$ = 39, 19.6 and 11.5 GeV pidv2_bes1 ;
pidprl_bes1 . The middle panel: The $v_{2}$ of quarks as a function of
$KE_{T}^{q}$= $\sqrt{(p_{T}/n_{q})^{2}+m_{q}^{2}}-m_{q}$, where $p_{T}$ is the
transverse momentum of hadron. Here we consider $m_{q}$ = 4 MeV/c2 and 100
MeV/c2 for light and strange quarks, respectively. The right panel: same as
the middle panel , but $m_{q}$ = 140 MeV/c2 for strange quarks. Blue dashed
curves are the fit to light quarks $v_{2}$ using Eq. 5 and red dashed lines
are the constant polynomial fit to the ratios shown in the respective bottom
panels.

###  II.2 The test of scaling at LHC energy

We have tested the scaling at LHC energy, where we have used the latest
measurement of $\phi$ and proton $v_{2}$ at $\sqrt{s_{NN}}$ = 5.02 TeV by
ALICE experiments alice_502 . We do not have data for $\Omega$ and $\Xi$ at
this energies. Figure 5 shows $v_{2}/n_{q}$ as a function of $p_{T}/n_{q}$ and
$KE_{T}^{q}$ in 20-30% Pb+Pb collisions at $\sqrt{s_{NN}}$= 5.02 TeV. We can
see that $v_{2}/n_{q}$ for $\phi$ meson is lower than that of protons at the
low $p_{T}$ region when plotted as a function of $p_{T}/n_{q}$. However, the
$v_{2}/n_{q}$ as a function of $KE_{T}^{q}$ shows similar values of $v_{2}$
for both proton and $\phi$. This is true for other centrality collisions as
well (not shown here). It is important to mention that the empirical scaling
of $v_{2}/n_{q}$ as a function of $KE_{T}^{q}$ works well at both RHIC and LHC
energies, where center-of-mass energies are very different, LHC is a factor
100 times higher than RHIC BES energies. This observation is interesting and
needs further theoretical studies to understand the true mechanism behind the
observed scaling in the measured hadrons $v_{2}$.

###  II.3 The test of scaling in higher order flow harmonics

In the end, we performed another systematic check using triangular flow
harmonic ($v_{3}$). It was shown before that, in case of $v_{3}$, one needs to
divide the $v_{3}$ by $n_{q}^{3/2}$ to see the scaling XuSun . Figure 6 (a)
shows $v_{3}/n_{q}^{3/2}$ as a function of $p_{T}/n_{q}$ for proton and $\phi$
measured by the STAR experiment in Au+Au collisions at $\sqrt{s_{NN}}$ = 200
GeV v3_200 . Like $v_{2}$, we do see the scaled $v_{3}$ of $\phi$ is lower
than scaled proton $v_{2}$ at low $p_{T}$. However, if we plot
$v_{3}/n_{q}^{3/2}$ as a function of $KE_{T}^{q}$, we do see good consistency
as observed for $v_{2}$.

Figure 5: (Color online) Panel (a) : The $v_{2}$ of light and strange quarks
obtained from proton, and $\phi$ as a function of $p_{T}/n_{q}$ in 20-30%
central Pb+Pb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV alice_502 . Panel (b) :
The $v_{2}$ of quarks as a function of $KE_{T}^{q}$=
$\sqrt{(p_{T}/n_{q})^{2}+m_{q}^{2}}-m_{q}$, where $p_{T}$ is the transverse
momentum of hadron. Here we consider $m_{q}$ = 4 MeV/c2 and 100 MeV/c2 for
light and strange quarks, respectively. Panel (c) : Same as panel (b), but
here $m_{q}$ = 140 MeV/c2 for strange quarks. Blue dashed curves are the fit
to light quarks $v_{2}$ using Eq. 5 and red dashed lines are the constant
polynomial fit to the ratios shown in the respective bottom panels.  Figure 6:
(Color online) Panel (a): The $v_{3}/n_{q}^{3/2}$ of proton and $\phi$ as a
function of $p_{T}/n_{q}$ in Au+Au collisions at $\sqrt{s_{NN}}$ = 200 GeV
v3_200 . Panel (b): $v_{3}/n_{q}^{3/2}$ of proton and $\phi$ a function of
$KE_{T}^{q}$= $\sqrt{(p_{T}/n_{q})^{2}+m_{q}^{2}}-m_{q}$, where $p_{T}$ is the
transverse momentum of hadron. Here we consider $m_{q}$ = 4 MeV/c2 and 100
MeV/c2 for light and strange quarks, respectively. Panel (c) : Same as panel
(b), but here $m_{q}$ = 140 MeV/c2 for strange quarks. Blue dashed curves are
the fit to light quarks $v_{3}$ using Eq. 5 and red dashed lines are the
constant polynomial fit to the ratios shown in the respective bottom panels.

##  III Summary

In summary, a compilation of the available data at RHIC and LHC energies for
the elliptic flow and triangular flow of identified hadrons is presented.
Within the framework of the quark recombination model, $v_{2}$ of light and
strange quarks are obtained from the measured $v_{2}$ of identified hadrons.
We find that the $v_{2}$ of strange quarks obtained from $\phi$-meson $v_{2}$
is found to be similar to that of from $\Omega$ $v_{2}$ in Au+Au collisions at
$\sqrt{s_{NN}}$ = 200 GeV, indicating that both $\phi$ and $\Omega$ are
produced through quark recombination at top RHIC energy. However, we find the
$v_{2}$ of strange quarks is consistently smaller than that of light quarks at
the low transverse momentum. We have shown an empirical relation between
$v_{2}$ of light and strange quarks, where $v_{2}$ of quarks scales as a
function of the transverse kinetic energy of quarks. The transverse kinetic
energy of quarks has been calculated using the bare mass of quarks. Such
empirical relation holds for both RHIC and LHC energies. It is shown that
$v_{3}$ of the measured hadrons also follows the same scaling. The observation
of such new scaling in flow measurement is purely empirical and needs further
theoretical investigation.

##  IV References

## References

  * (1) STAR Collaboration, J. Adams et al., Nucl. Phys. A 757, 102 (2005).
  * (2) J. Y. Ollitrault, Phys. Rev. D 46, 229 (1992); H. Sorge, Phys. Rev. Lett. 78, 2309 (1997).
  * (3) P. Huovinen et al., Phys. Lett. B 503, 58 (2001).
  * (4) P. F. Kolb et al., Phys. Rev. C 62, 054909 (2000).
  * (5) S. A. Voloshin et al., arXiv:0809.2949.
  * (6) STAR Collaboration, L. Adamczyk et al., Phys. Rev. C 88, 014902 (2013).
  * (7) STAR Collaboration, L. Adamczyk et al., Phys. Rev. Lett. 110, 142301 (2013).
  * (8) STAR Collaboration, L. Adamczyk et al., Phys. Rev. Lett. 116, 062301 (2016).
  * (9) STAR Collaboration, J. Adams et al., Phys. Rev. Lett. 95, 122301 (2005); PHENIX Collaboration, A. Adare et al., Phys. Rev. Lett. 98, 162301 (2007).
  * (10) PHOBOS Collaboration, B. B. Back et al., Phys. Rev. Lett. 94, 122303 (2005).
  * (11) STAR Collaboration, J. Adams et al., Phys. Rev. C 72, 014904 (2005); B. I. Abelev et al., Phys. Rev. C 75, 054906 (2007); B. I. Abelev et al., Phys. Rev. C 77, 054901 (2007).
  * (12) PHENIX Collaboration, S. Afanasiev et al., Phys. Rev. C 80, 024909 (2009); A. Adare et al., Phys. Rev. Lett. 105, 062301 (2010); S. Afanasiev et al., Phys. Rev. Lett. 99, 052301 (2007).
  * (13) R.J.Fries et al., Phys.Rev.C 68, 044902 (2003).
  * (14) P. F. Kolb et al., Phys.Rev. C 69, 051901 (2004).
  * (15) D. Molnar, S. A. Voloshin; Phys. Rev. Lett. 91, 092301 (2003).
  * (16) L.-W. Chen and C. M. Ko. Phys. Rev. C 73, 044903 (2006).
  * (17) J. H. Chen et al., Phys. Rev. C 74, 064902 (2006).
  * (18) S. Choudhury, et al., Phys. Rev. C 95, 024904, (2017).
  * (19) J. Song, et al., Phys. Rev. C 103, 034907 (2021).
  * (20) R.J.Fries et al., Ann. Rev. Nucl. Part. Sci. 58, 177-205, (2008).
  * (21) STAR Collaboration, L. Adamczyk et al., Phys. Rev. C 93, 021903 (2016).
  * (22) J. Song, et al., Eur. Phys. J. C 81, 5 (2021).
  * (23) R.Fries, et al., Ann. Rev. Nucl. Part. Sci. 58:177-205 (2008).
  * (24) J. Jia, C.Zhang Phys. Rev. C 75: 031901 (2007).
  * (25) PHENIX Collaboration, A. Adare et al., Phys. Rev. Lett.98,162301,(2007).
  * (26) P. A. Zyla et al., (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020) and 2021 update.
  * (27) ALICE Collaboration, S. Acharya et al., JHEP 09, 006 (2018).
  * (28) X. Sun (for the STAR Collaboration) Nucl. Phys. A 931, 1194-1198 (2014).
  * (29) STAR Collaboration, M. S. Abdallah et al., arXiv:2203.07204.