\section{Introduction}

Algorithm designs for combinatorial optimization problems (COPs) are important and challenging tasks. A wide variety of real-world problems are COPs, such as vehicle routing \cite{toth2002vehicle}, network design \cite{johnson1978complexity}, path planning \cite{pohl1970heuristic} and mechanism design \cite{de2003combinatorial} problems, and a majority of them are NP-hard to solve. In the past few decades, algorithms, including optimal algorithms, approximation algorithms and heuristic algorithms, have been studied extensively due to the importance of COPs. Those algorithms are mostly designed by human through costly processes that often require deep understanding of the problem domains and their underlying structures as well as considerable time and effort.

Recently, there has been an increased interest to automate algorithm designs for COPs with machine learning (ML). Many  ML approaches learn to either construct  or improve solutions within an algorithmic framework, such as greedy search, local search or tree search, for a specific COP, such as the traveling salesman problem (TSP) \cite{xin2021neurolkh,zheng2021combining}, vehicle routing problem
 (VRP) \cite{kool2018attention} or independent set problem \cite{li2018combinatorial}, and are often not easily applicable to other COPs.

In contrast, Integer Linear Programs (ILPs) can flexibly encode and solve a broad family of COPs, such as minimum vertex cover, set covering and facility location problems. ILPs can be solved by Branch and Bound (BnB) \cite{land2010automatic}, an optimal tree search algorithm that can achieve state-of-the-art for ILPs. Over the past decades, BnB has been improved tremendously to become the core of many popular ILP solvers such as SCIP \cite{BestuzhevaEtal2021OO}, CPLEX \cite{cplex2009v12} and Gurobi \cite{gurobi}. However, due to its exhaustive search nature, it is hard for BnB to scale to large instances \cite{khalil2016learning,gasse2019exact}.

On the other hand, Large Neighborhood Search (LNS) has recently been shown to  find high quality solutions much faster than BnB for large ILP instances \cite{song2020general,wu2021learning,sonnerat2021learning,huang2022local}. LNS starts from an initial solution (i.e., a feasible assignment of values to variables) and then improves the current best solution by iteratively picking a subset of variables to reoptimize while leaving others fixed. Picking which subset to reoptimize, i.e., the \emph{destroy heuristic}, is a critical component in LNS. Hand-crafted destroy heuristics, such as the randomized heuristic \cite{song2020general,sonnerat2021learning} and the Local Branching (LB) heuristic \cite{fischetti2003local}, are often either inefficient (slow to find good subsets) or ineffective (find subsets of bad quality). ML-based destroy heuristics have also been proposed and outperform hand-crafted ones. State-of-the-art approaches include \DM~\citep{sonnerat2021learning} that uses imitation learning (IL) to imitates the LB heuristic and \RL~\citep{wu2021learning} that uses a similar framework to \DM but trained with reinforcement learning (RL).

In this paper, we propose a novel ML-based LNS for ILPs, namely \emph{\CL}, that uses contrastive learning (CL) \cite{chen2020simple,khosla2020supervised} to learn efficient and effective destroy heuristics. Similar to \DM~\cite{sonnerat2021learning}, we learn to imitate the \emph{Local Branching (LB)} heuristic, a destroy heuristic that selects the optimal subset of variables within the Hamming ball of the incumbent solutions. LB requires
solving another ILP with same size as the original problem and thus is computationally expensive. We not only use the optimal subsets provided by LB as the expert demonstration (as in \DM), but also leverage intermediate solutions and perturbations. When solving the ILP for LB, intermediate solutions are found and those that are close to optimal in term of effectiveness become \emph{positive samples}.
We also collect \emph{negative samples} by randomly perturbing the optimal subset.
With both positive and negative samples, instead of a classification loss as in \DM, we use a contrastive loss that encourages the model to predict the subset similar to the positive samples but dissimilar to the negative ones with similarity measured by dot products \cite{oord2018representation,he2020momentum}.
Finally, we also use a richer set of features and use graph attention networks (GAT) instead of GCN, to further boost performance.

Empirically, we show that \CL outperforms state-of-the-art ML and non-ML approaches at different runtime cutoffs ranging from a few minutes to an hour in terms of multiple metrics, including the primal gap, the primal integral, the best performing rate and the survival rate, demonstrating the effectiveness and efficiency of \CL. In addition, \CL shows great generalization performance on test instances two times larger than training instances.

\section{Background}
In this section, we first define ILPs and then introduce LNS for ILP solving and the Local Branching (LB) heuristic.

\subsection{ILPs}

An \textit{integer linear program (ILP)} is defined as
\begin{equation}
    \min \bc^{\mathsf{T}}\bx \textrm{ s.t. } \bA \bx\leq \bb \textrm{ and } \bx \in \{0,1\}^n, \label{eqn::ILPformulation}
\end{equation}
where $\bx = (x_1,\ldots, x_n)^\sfT$ denotes the $n$ binary variables to be optimized, $\bc\in \mathbb{R}^n$ is the vector of objective coefficients, $\bA\in \mathbb{R}^{m\times n}$ and $\bb\in \mathbb{R}^{m}$ specify $m$ linear constraints.
A \textit{solution} to the ILP is a feasible assignment of values to the variables.
In this paper, we focus on the formulation above that consists of only binary variables, but our methods can be applied to mixed integer linear programs with continuous variables and/or non-binary integer variables.

\iffalse
\begin{algorithm}[t]
\small
	\caption{LNS for ILPs}\label{algo::LNSforILP}
	\begin{algorithmic}[1]
	     \State {\bf Input: } An ILP.
	     \State $\bx^0\gets$ Find an intial solution to the input ILP
	     \State $t\gets 0$
	    \While {time limit not exceeded}
		    \State $\calX^t\gets$ Select a subset of variables to destroy
		    \State $\bx^{t+1}\gets$ Solve the ILP with additional constraints $\{x_i = x_i^t: x_i\notin \calX^t\}$
		    \State $t\gets t + 1$

		\EndWhile
		\State \Return $\bx^t$
	\end{algorithmic}
\end{algorithm}
\fi

\subsection{LNS for ILP solving}\label{sec::LNSforILPs}

LNS is a heuristic algorithm that starts with an initial solution and then iteratively destroys and reoptimizes a part of the solution until a runtime limit is exceeded or some stopping condition is met. Let $\mathcal{I}=(\bA, \bb, \bc)$ be the input ILP, where $\bA, \bb$ and  $\bc$ are the coefficients defined in Equation (\ref{eqn::ILPformulation}), and $\bx^0$ be the initial solution (typically found by running BnB for a short runtime). In iteration $t\geq 0$ of LNS, given the \textit{incumbent solution} $\bx^t$,  defined as the best solution found so far, a \textit{destroy heuristic} selects a subset of $k^t$ variables $\calX^t= \{x_{i_1},\ldots, x_{i_{k^t}}\}$. The reoptimization is done by solving a sub-ILP with $\calX^t$ being the variables while fixing the values of $x_j\notin \calX^t$ the same as in $\bx^t$. The solution to the sub-ILP is the new incumbent solution $\bx^{t+1}$ and then LNS proceeds to iteration $t+1$.
Compared to BnB, LNS is more effective in improving the objective value $\bc^{\mathsf{T}}x$ especially on difficult instances \cite{song2020general,sonnerat2021learning,wu2021learning}. Compared to other local search methods, LNS explores a large neighborhood in each step and thus, is more effective in avoiding local minima.

 \paragraph{Adaptive Neighborhood Size}
Adaptive methods are commonly used to set the neighborhood size $k^t$ in previous work \cite{sonnerat2021learning,huang2022local}. The initial neighborhood size $k^0$ is set to a constant or a fraction of the number of variables. In this paper, we consider the following adaptive method \cite{huang2022local}: in iteration $t$, if LNS finds an improved solution, we let $k^{t+1}=k^t$, otherwise $k^{t+1} =\min\{ \gamma\cdot k^t, \beta\cdot n\}$ where  $\gamma>1$ is a constant and we upper bound  $k^t$ to a constant fraction $\beta<1$ of the number of variables to make sure the sub-ILP is not too large (thus, too difficult) to solve. Adaptively setting $k^t$ helps LNS escape local minima by expanding the search neighborhood when it fails to improve the solution.

\subsection{LB Heuristic}

The LB Heuristic \cite{fischetti2003local} is originally proposed as a primal heuristic in BnB but also applicable in LNS for ILP solving \cite{sonnerat2021learning,liu2022learning}.
Given the incumbent solution $\bx^t$ in iteration $t$ of LNS, LB aims to find the subset of variables to destroy $\calX^t$ such that it leads to the optimal $\bx^{t+1}$ that differs from $\bx^t$ on at most $k^t$ variables, i.e., it computes the optimal solution $\bx^{t+1}$ that sits within a given Hamming ball of radius $k^t$ centered around $\bx^t$.
To find $\bx^{t+1}$, the LB heuristic solves the LB ILP that is exactly the same ILP from input but with one additional constraint that limits the distance between $\bx^t$ and $\bx^{t+1}$:
$\sum_{i\in[n]:x^t_i=0}x^{t+1}_i + \sum_{i\in[n]:x^t_i=1}(1-x^{t+1}_i)\leq k^t. $
The LB ILP is of the same size of the input ILP (i.e., it has the same number of variables and one more constraint), therefore, it is often too slow to be useful in practice.

\section{Related Work}

In this section, we summarize related work on LNS for ILPs and other COPs, learning to solve ILPs with BnB and contrastive learning for COPs. We also summarize additional related work on LNS-based primal heuristics for BnB and learning to solve other COPs in Appendix.

\subsection{LNS for ILPs and Other COPs}
Huge effort has been made to improve BnB for ILPs in the past decades, but LNS for ILPs has not been studied extensively.

Recently, \citet{song2020general} show that even a randomized destroy heuristic in LNS can outperform state-of-the-art BnB. They also show that an ML-guided decomposition-based LNS can achieve even better performance, where they apply RL and IL to learn destroy heuristics that decompose the set of variables into equally-sized subsets using a classification loss.
\citet{sonnerat2021learning} learn to select variables by imitating LB. \RL~\citep{wu2021learning} uses a similar framework but trained with RL and outperforms \citet{song2020general}.
Both \citet{wu2021learning} and \citet{sonnerat2021learning} use the bipartite graph representations of ILPs to learn the destroy heuristics represented by GCNs.
Another line of related work focuses on improving LB. \citet{liu2022learning} use ML to tune the runtime limit and neighborhood sizes for LB. \citet{huang2022local} propose \LBRELAX to select variables by solving  the LP relaxation of LB.

Besides ILPs, LNS has been applied to solve many COPs, such as VRP \cite{ropke2006adaptive,azi2014adaptive}, TSP \cite{smith2017glns}, scheduling \cite{kovacs2012adaptive,vzulj2018hybrid} and path planning problems \cite{LiAAAI22,LiIJCAI21}. ML methods have also been applied to improve LNS for those applications \cite{chen2019learning,lu2019learning,hottung2020neural,li2021learning,huang2022anytime}.

\subsection{Learning to Solve ILPs with BnB}

Several studies have applied ML to improve BnB. The majority of works focus on learning to either  select variables to branch on \cite{khalil2016learning,gasse2019exact,gupta2020hybrid,zarpellon2021parameterizing} or select nodes to expand \cite{he2014learning,labassi2022learning}. There are also works on learning to schedule and run primal heuristics \cite{khalil2017learning,chmiela2021learning} and to select cutting planes \cite{tang2020reinforcement,paulus2022learning,huang2022learning}.

\subsection{Contrastive Learning for COPs}

While contrastive learning of visual representations \cite{hjelm2018learning,he2020momentum,chen2020simple} and graph representations \cite{you2020graph,tong2021directed}
have been studied extensively, it has not been explored much for COPs. \citet{mulamba2021contrastive} derive a contrastive loss for decision-focused learning to solve COPs with uncertain inputs that can be learned from historical data, where they view non-optimal solutions as negative samples. \citet{duan2022augment} use contrastive pre-training to learn good representations for the boolean satisfiability problem.

\section{Contrastive Learning for LNS}

Our goal is to learn a policy, a destroy heuristic represented by an ML model, that selects a subset of variables to destroy and reoptimize in each LNS iteration. Specifically, let $\bs^t=(\calI, \bx^t)$ be the current state in iteration $t$ of LNS where $\mathcal{I}=(\bA, \bb, \bc)$ is  the ILP and $\bx^t$ is the incumbent solution, the policy predicts an action $\ba^t=(a_1^t,\ldots, a_n^t)\in \{0,1\}^n$, a binary representation of the selected variables $\calX^t$
indicating whether $x_i$ is selected ($a^t_i=1$) or not ($a^t_i=0$). We use contrastive learning to learn to predict high quality $\ba^t$ such that, after solving the sub-ILP derived from $\ba^t$ (or $\calX^t$), the resulting incumbent solution $\bx^{t+1}$ is improved as much as possible.
Next, we describe how we prepare data for contrastive learning, the policy network and the contrastive loss used in training, and finally introduce how the learned policy is used in \CL.

\iffalse
\subsection{MDP Formulation}

 Our goal is to learn a policy, a destroy heuristic represented by an ML model, that selects a subset of variables to destroy and reoptimize in each LNS iteration. We formulate this sequential decision problem as an MDP where step $t\geq 0$ of the MDP corresponds to iteration $t$ in LNS.
The state $\bs^t = (\mathcal{I},\bx^t)$ at step $t\geq0$ of the MDP consists of the ILP $\mathcal{I}=(\bA, \bb, \bc)$ and its incumbent solution $\bx^t$. The action $\ba^t=(a_1^t,\ldots, a_n^t)\in \{0,1\}^n$ at step $t$ is a binary representation of the selected variables $\calX^t$ where $a^t_i = \mathbf{1}_{x_i\in \calX^t}$ indicates whether $x_i$ is selected or not.
The transition from  $\bs^t$ with action $\ba^t$ to $\bs^{t+1}=(\mathcal{I},\bx^{t+1})$ is done by solving the sub-ILP that reoptimizes $\calX^t$ to obtain $\bx^{t+1}$.
The reward $r^t=r(\bs^t,\ba^t)= \bc^{\mathsf{T}}(\bx^{t}-\bx^{t+1})\geq 0$ which is the improvement of the objective value.

An episode begins with $\bs^0=(\mathcal{I},\bx^0)$, the input ILP and its initial solution, and proceeds by running the MDP until the runtime limit is exceeded or some stopping condition is met. Let $T$ be the index of the final MDP step, the goal is to maximize the total reward $\sum_{t\leq T} r^t$, which is equivalent to minimizing the objective value of the ILP.
\fi

\subsection{Data Collection}

Following previous work by \citet{sonnerat2021learning}, we use LB as the expert policy to collect good demonstrations to learn to imitate. Formally, for a given state $\bs^t=(\calI,\bx^t)$, we use LB to find the optimal action $\ba^t$ that leads to the minimum $\bc^{\mathsf{T}}\bx^{t+1}$ after solving the sub-ILP.
Different from the previous work, we use contrastive learning to learn to make discriminative predictions of  $\ba^t$ by contrasting positive and negative samples (i.e., good and bad examples of actions $\ba^t$). In the following, we describe how we collect the positive sample set $\calS_{\sfp}^t$ and the negative sample set $\calS_{\sfn}^t$.

\paragraph{Collecting Positive Samples $\calS_{\sfp}^t$}

During data collection, given $\bs^t=(\calI, \bx^t)$, we solve the LB ILP with the incumbent solution $\bx^t$ and neighborhood size $k^t$ to find the optimal $\bx^{t+1}$.
LNS proceeds to  iteration $t+1$ with $\bx^{t+1}$ until no improving solution $\bx^{t+1}$ could be found by the LB ILP within a runtime limit.
In experiments, the LB ILP is solved with SCIP 8.0.1 \cite{BestuzhevaEtal2021OO} with an hour runtime limit  and $k^t$ is fine-tuned for each type of instances. After each solve of the LB ILP, in addition to the best solution found, SCIP records all intermediate solutions found during the solve. We look for intermediate solutions $\bx'$ whose resulting improvements on the objective value is at least $0<\alpha_{\sfp}\leq 1$ times the best improvement  (i.e., $\bc^{\mathsf{T}}(\bx^t-\bx')\geq \alpha_{\sfp} \cdot \bc^{\mathsf{T}}(\bx^t-\bx^{t+1})$) and consider their corresponding actions as positive samples. We limit the number of the positive samples $|\calS^t_{\sfp}|$ to  $u_{\sfp}$.
If more than $u_{\sfp}$ positive samples are available, we record the top $u_{\sfp}$ ones to avoid large computational overhead with too many samples when computing the contrastive loss (see Section \ref{sec::trainingwithCL}).
$\alpha_{\sfp}$ and $u_{\sfp}$ are set to $0.5$ and $10$, respectively, in experiments.

\paragraph{Collecting Negative Samples $\calS_{\sfn}^t$}

Negative samples are critical parts of contrastive learning to help distinguish between good and bad demonstrations. We collect a set of $c_{\sfn}^t$ negative samples $\calS_{\sfn}^t$, where $c_{\sfn}^t=\kappa |\calS_{\sfp}^t|$ and $\kappa$ is a hyperparameter to control the ratio between the numbers of positive and negative samples.
Suppose $\calX^t$ is the optimal set of variables selected by LB. We then perturb $\calX^t$ to get $\hat{\calX}^{t}$ by replacing $5\%$ of the variables in  $\calX^t$ with the same number of those not in $\calX^t$ uniformly at random. We then solve the corresponding sub-ILP derived from $\hat{\calX}^{t}$ to get a new incumbent solution $\hat{\bx}^{t+1}$. If the resulting improvement of $\hat{\bx}^{t+1}$ is less than $0\leq \alpha_{\sfn}<1$ times the best improvement  (i.e., $\bc^{\mathsf{T}}(\bx^t-\hat{\bx}^{t+1})\leq \alpha_{\sfn} \cdot \bc^{\mathsf{T}}(\bx^t-\bx^{t+1})$), we consider its corresponding action as a negative sample. We repeat this $c_{\sfn}^t$ times to collect negative samples. If less than $c_{\sfn}^t$ negative samples is collected, we increase the perturbation rate from $5\%$ to $10\%$ and generate another $c_{\sfn}^t$ samples. We keep increasing the perturbation rate at an increment of $5\%$ until $c_{\sfn}^t$ negative samples are found or it reaches $100\%$. In experiments, we set $\kappa = 9$ and $\alpha_{\sfn}= 0.05$, and it takes less than 3 minutes to collect negative samples for each state.

\subsection{Policy Network}
\label{sec::policyNetwork}
Following previous work on learning for ILPs \cite{gasse2019exact, sonnerat2021learning, wu2021learning}, we use a bipartite graph representation of ILP to encode a state $\bs^{t}$. The bipartite graph consists of $n+m$ nodes representing the $n$ variables and $m$ constraints on two sides, respectively, with an edge connecting a variable and a constraint if the variable has a non-zero coefficient in the constraint. Following \citet{sonnerat2021learning}, we use features proposed in \citet{gasse2019exact} for node features and edge features in the bipartite graph and also include a fixed-size window of most recent incumbent values as variable node features with the window size set to 3 in experiments. In addition to features used in \citet{sonnerat2021learning}, we include features proposed in \citet{khalil2016learning} computed at the root node of BnB to make it a richer set of variable node features.

We learn a policy $\pi_{\btheta}(\cdot)$ represented by a graph attention network (GAT) \cite{brody2021attentive} parameterized by learnable weights $\btheta$. The policy takes as input the state $\bs^t$ and output a score vector $\pi_{\btheta}(\bs^t)\in [0,1]^n$, one score per variable.
To increase the modeling capacity and to manipulate node interactions proposed by our architecture, we use embedding layers to map each node feature and edge feature to space $\mathbb{R}^d$.
Let $\bfv_j, \bfc_i, \bfe_{i,j}\in\mathbb{R}^d$ be the embeddings of the $i$-th variable, $j$-th constraint and the edge connecting them output by the embedding layers. Since our graph is bipartite, following previous work \cite{gasse2019exact}, we perform two rounds of message passing through the GAT.  In the first round, each constraint node $\bfc_i$ attends to its neighbors $\calN_i$ using an attention structure with $H$ attention heads
to get updated constraint embeddings $\bfc'_i$ (computed as a function of $\bfv_j, \bfc_i, \bfe_{i,j}$).
In the second round, similarly, each variable node attends to its neighbors to get  updated variable embeddings $\bfv'$ (computed as a function of $\bfv_j, \bfc'_i, \bfe_{i,j}$) with another set of attention weights.
After the two rounds of message passing, the final representations of variables $\bfv'$ are passed through a multi-layer perceptron (MLP) to obtain a scalar value for each variable and, finally, we apply the sigmoid  function to get a score  between 0 and 1.
Full details of the network architecture are provided in Appendix.
In experiments, $d$ and $H$ are set to $64$ and $8$, respectively.

\begin{table*}[htbp]
\centering
\caption{ Names and the average numbers of variables and constraints of the test instances.
\label{table::instanceSize}}
\scriptsize
\begin{tabular}{c|crrr|crrr}
\hline
              & \multicolumn{4}{c|}{Small Instances}                                                                              & \multicolumn{4}{c}{Large Instances}                                                                                \\ \hline
Name          & \multicolumn{1}{c}{MVC-S}  & \multicolumn{1}{c}{MIS-S} & \multicolumn{1}{c}{CA-S}  & \multicolumn{1}{c|}{SC-S} & \multicolumn{1}{c}{MVC-L}   & \multicolumn{1}{c}{MIS-L}  & \multicolumn{1}{c}{CA-L}  & \multicolumn{1}{c}{SC-L} \\ \hline
\#Variables   & \multicolumn{1}{r}{1,000}  & \multicolumn{1}{r}{6,000} & \multicolumn{1}{r}{4,000} & 4,000                     & \multicolumn{1}{r}{2,000}   & \multicolumn{1}{r}{12,000} & \multicolumn{1}{r}{8,000} & 8,000                    \\
\#Constraints & \multicolumn{1}{r}{65,100} & \multicolumn{1}{r}{23,977}      & \multicolumn{1}{r}{2,675}      & 5,000                     & \multicolumn{1}{r}{135,100} & \multicolumn{1}{r}{48,027}       & \multicolumn{1}{r}{5,353}      & 5,000                    \\ \hline
\end{tabular}
\end{table*}

\subsection{Training with a Contrastive Loss}\label{sec::trainingwithCL}

Given a set of ILP instance for training, we follow the expert's trajectory to collect training data. Let $\calD = \{(\bs, \calS_{\sfp}, \calS_{\sfn})\}$ be the set of states with their corresponding sets of positive and negative samples in the training data.
A contrastive loss is a function whose value is low when the predicted action $\pi_{\btheta}(\bs)$ is similar to the positive samples $\calS_{\sfp}$ and dissimilar to the negative samples $\calS_{\sfn}$. With similarity measured by dot products, a form of supervised contrastive loss, called InfoNCE \cite{oord2018representation,he2020momentum},
is used in this paper:
$$\small
\calL(\btheta) = \sum_{(\bs, \calS_{\sfp}, \calS_{\sfn})\in \calD}\frac{-1}{|\calS_{\sfp}|}\sum_{\ba\in \calS_{\sfp}}\log\frac{\exp(\ba^{\mathsf{T}}\pi_{\btheta}(\bs)/\tau)}{\sum_{\ba'\in \calS_{\sfn}\cup{\{\ba\}}}\exp(\ba'^{\mathsf{T}}\pi_{\btheta}(\bs)/\tau)}
$$
where  $\tau$ is a temperature hyperparameter set to 0.07 \cite{he2020momentum} in experiments.

\begin{figure*}[tbp]
    \centering
    \includegraphics[width=0.8\textwidth]{figure_all/legend_ML_horizontal_timeVSobj_3601.eps}

    \subfloat[MVC-S (left) and MVC-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSgap_MVC_BA70_1000_3601_nh_50.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSgap_MVC_BA70_2000_3601_nh_50.eps}
    }
    \subfloat[MIS-S (left) and MIS-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSgap_INDSET_ER4_6000_3601_nh_3000.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSgap_INDSET_ER4_12000_3601_nh_6000.eps}
    }\\
    \subfloat[CA-S (left) and CA-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSgap_CAT_2000_4000_0.7_3601_nh_1000.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSgap_CAT_4000_8000_0.7_3601_nh_1000.eps}
    }
    \subfloat[SC-S (left) and SC-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSgap_SC_4000_5000_3601_nh_150.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSgap_SC_8000_5000_3601_nh_150.eps}
    }
    \caption{The primal gap (the lower the better) as a function of runtime, averaged over 100 test instances. For ML approaches, the policies are trained on only small training instances but tested on both small and large test instances. \label{res::gap}}
\end{figure*}

\subsection{Applying Learned Policy $\pi_{\btheta}$} \label{sec::usingthepolicy}
We apply the learned policy $\pi_{\btheta}$ in LNS.
In iteration $t$, let $(v_1,\cdots,v_n):=\pi_{\btheta}(\bs^t)$ be the variable scores output by the policy. To select $k^t$ variables, \CL greedily selects those with the highest scores. Previous works \cite{sonnerat2021learning,wu2021learning} commonly use sampling methods to select the variables, but those sampling methods are empirically worse than our greedy method in \CL. However, when the adaptive neighborhood size $k^t$ reaches its upper bound $\beta\cdot n$, \CL may repeat the same prediction due to deterministic selection process. When this happens, we switch to the sampling method introduced in \cite{sonnerat2021learning}. The sampling method selects variables sequentially: at each step, a variable $x_i$ that has not been selected yet is selected with probability proportional to $v_i^{\eta}$, where $\eta$ is a temperature parameter set to $0.5$ in experiments.

\section{Empirical Evaluation}

In this section, we introduce our evaluation setup and then present the results. Our code will be made available to the public upon publication.

\subsection{Setup}

\paragraph{Instance Generation} We evaluate on four NP-hard problem benchmarks that are widely used in existing studies \cite{wu2021learning,song2020general,scavuzzo2022learning}, which consist of two graph optimization problems, namely the minimum vertex cover (MVC) and maximum independent set (MIS) problems, and two non-graph optimization problems, namely the combinatorial auction (CA) and set covering (SC) problems.
We first generate a test set of 100 \emph{small instances} for each problem, namely MVC-S, MIS-S, CA-S and SC-S.
MVC-S instances are generated according to the Barabasi-Albert random graph model \cite{albert2002statistical}, with 1,000 nodes and average degree 70 following \cite{song2020general}.
MIS-S instances are generated according to the Erdos-Renyi random graph model \cite{erdos1960evolution}, with 6,000 nodes and average degree 5 following \cite{song2020general}.
CA-S instances are generated with 2,000 items and 4,000 bids according to the arbitrary relations in \citet{leyton2000towards}.
SC-S instances are generated with 4,000 variables and 5,000 constraints following \citet{wu2021learning}.
We then generate another test set of 100 \emph{large instances} for each problem by doubling the number of variables, namely MVC-L, MIS-L, CA-L and SC-L.
For each test set,  Table \ref{table::instanceSize} shows its average numbers of variables and constraints. More details of instance generation are included in Appendix.

 For data collection and training, we generate another set of 1,024 small instances for each problem. We split these instances into training and validation sets, each consisting of 896 and 128 instances, respectively.

\paragraph{Baselines}

We compare \CL with five baselines:
    (1) BnB: using SCIP (v8.0.1), the state-of-the-art open-source ILP solver, with the aggressive mode fine-tuned to focus on improving the objective value;

    (2) \RANDOM: LNS which selects the neighborhood by uniformly sampling $k^t$ variables without replacement;

    (3) \LBRELAX \cite{huang2022local}: LNS which selects the neighborhood with the \LBRELAX heuristics;
    (4) \DM \cite{sonnerat2021learning};
    (5) \RL \cite{wu2021learning}. We compare with two more baselines in Appendix.
\iffalse
\begin{itemize}
    \item BnB using SCIP (v8.0.1), the state-of-the-art open-source ILP solver, with the aggressive mode turned on to focus on improving the objective value;

    \item \RANDOM: LNS which selects the neighborhood by uniformly sampling $k^t$ variables without replacement;

    \item \LBRELAX \cite{huang2022local}: LNS which selects the neighborhood with the \LBRELAX heuristics;

    \item \DM \cite{sonnerat2021learning}: We implement \DM since the authors do not fully open source the code and fine tune its hyperparameters for each problem;
    \item \RL \cite{wu2021learning}: This approach does not require a given neighborhood size $k^t$ since it is defined implicitly by how the policy is used. We use the code and hyperparameters made available by the authors.
\end{itemize}
\fi
For each ML approach,  a separate model is trained for each problem on the small training set and tested on both small and large test sets. We implement \DM  and fine tune its hyperparameters for each problem since the authors do not fully open source the code. \DM uses the same training dataset as \CL but uses only the positive samples. For \RL, we use the code and hyperparameters provided by the authors and train the models with five random seeds to select one with the best performance on the validation sets. We do not compare to the approach by \citet{song2020general} since it performs worse than \RL on multiple problems \cite{wu2021learning}.

\paragraph{Metrics}
We use the following metrics to evaluate all approaches:
(1) The \textit{primal bound} is the objective value of the ILP;
(2) The \textit{primal gap} \cite{berthold2006primal} is the normalized difference between the primal bound $v$ and a precomputed best known objective value $v^*$, defined as $\frac{|v-v^*|}{\max(v,v^*,\epsilon)}$ if $v$ exists and $v\cdot v^*\geq 0$, or 1 otherwise. We use $\epsilon=10^{-8}$ to avoid division by zero;
(3) The \textit{primal integral} \cite{achterberg2012rounding} at time $q$  is the integral on $[0,q]$ of the primal gap as a function of runtime. It captures the quality of and the speed at which solutions are found;
(4) The \textit{survival rate}  to meet a certain primal gap threshold is the fraction of instances with primal gaps below the threshold \cite{sonnerat2021learning}; (5) The \emph{best performing rate}  of an approach is  the fraction of instances on which it achieves the best primal gap (including ties) compared to all approaches at a given runtime cutoff.
Since BnB and LNS are both anytime algorithms, we  show these metrics as a function of runtime or the number of iterations in LNS (when applicable) to demonstrate their anytime performance.

\begin{table}[bp]
\scriptsize
\centering
\caption
Primal gap (PG) (in percent), primal integral (PI) at 60 minutes runtime cutoff, averaged over 100 test instances and their standard deviations. ``$\downarrow$'' means the lower the better. For ML approaches, the policies are trained on only small training instances but tested on both small and large test instances. \label{res::smalltable60}}

\begin{tabular}{c|rr|rr}
\hline
                                   & \multicolumn{1}{c}{PG (\%) $\downarrow$}         & \multicolumn{1}{c|}{PI $\downarrow$}               & \multicolumn{1}{c}{PG (\%) $\downarrow$}         & \multicolumn{1}{c}{PI $\downarrow$} \\ \hline

                       & \multicolumn{2}{c|}{MVC-S}                                                                        & \multicolumn{2}{c}{MIS-S}                                                                         \\ \hline
BnB                     & {1.32$\pm$0.43} & 66.1$\pm$13.1            & {5.10$\pm$0.69} & 222.8$\pm$25.9           \\
\RANDOM                & {0.96$\pm$1.26} & 38.0$\pm$44.8             & {0.24$\pm$0.14} & 22.1$\pm$5.0             \\
\LBRELAX                 & {1.38$\pm$1.51} & 57.0$\pm$51.2             & {0.65$\pm$0.20} & 46.9$\pm$6.5             \\
\DM        & {0.29$\pm$0.23} & {19.2$\pm$10.2}               & {0.22$\pm$0.17} & { 19.4$\pm$5.8}             \\
\RL               & {0.61$\pm$0.34} & 29.6$\pm$11.5             & {0.22$\pm$0.14} & 17.2$\pm$5.2            \\
{\bf\CL}                & {\bf 0.17$\pm$0.09} & {\bf 8.7$\pm$6.7}             & {\bf 0.15$\pm$0.15} & {\bf12.8$\pm$5.4}            \\ \hline
\multicolumn{1}{l|}{} & \multicolumn{2}{c|}{CA-S}                                                                         & \multicolumn{2}{c}{SC-S}                                                                          \\ \hline
BnB                & 2.28$\pm$0.59 & 137.4$\pm$25.9      & 1.13$\pm$0.95 & 86.7$\pm$37.9                         \\
\RANDOM             & 5.90$\pm$1.02 & 235.6$\pm$34.9      & 2.67$\pm$1.29 & 124.3$\pm$45.4                    \\
\LBRELAX       & 1.65$\pm$0.57 & 140.5$\pm$18.3      & 0.86$\pm$0.83 & 63.2$\pm$31.6    \\
\DM             &1.09$\pm$0.51 & 90.0$\pm$20.8   & 1.33$\pm$0.97 & 63.2$\pm$34.3            \\
\RL            & 6.32$\pm$1.03 & 249.2$\pm$35.9   & 1.10$\pm$0.77 & 77.8$\pm$28.9            \\
{\bf\CL} & {\bf 0.65$\pm$0.32} & {\bf50.7$\pm$22.7} & {\bf0.50$\pm$0.58} & {\bf26.2$\pm$12.8}\\ \hline
\hline

                       & \multicolumn{2}{c|}{MVC-L}                                                                        & \multicolumn{2}{c}{MIS-L}                                                                         \\ \hline
BnB                    & {2.41$\pm$0.40} & 130.2$\pm$11.1            & {6.29$\pm$1.62} & 285.1$\pm$18.2           \\
\RANDOM                 & {0.38$\pm$0.24} & 22.7$\pm$8.0              & {\bf0.11$\pm$0.08} & 19.0$\pm$3.1             \\
\LBRELAX                 & {0.46$\pm$0.23} & 48.4$\pm$7.5             & {0.91$\pm$0.16} & 68.6$\pm$5.5             \\
\DM               & {0.27$\pm$0.23} & {21.2$\pm$8.1}               & {0.29$\pm$0.15} & { 27.1$\pm$5.5}             \\
\RL               & {0.59$\pm$0.30} & 37.3$\pm$9.6             & {0.14$\pm$0.12} & 18.9$\pm$4.1            \\
{\bf\CL}                & {\bf 0.05$\pm$0.04} & {\bf 9.1$\pm$3.4}             & { 0.12$\pm$0.11} & {\bf12.9$\pm$4.4}            \\ \hline
\multicolumn{1}{l|}{} & \multicolumn{2}{c|}{CA-L}                                                                         & \multicolumn{2}{c}{SC-L}                                                                          \\ \hline
BnB        & 2.74$\pm$1.87 & 320.9$\pm$83.1           & 1.54$\pm$1.33 & 115.0$\pm$42.5            \\
\RANDOM     & 5.37$\pm$0.75 & 229.2$\pm$24.4          & 3.31$\pm$1.79 & 166.4$\pm$61.3         \\
\LBRELAX    & 1.61$\pm$1.50 & 153.0$\pm$50.3          & 1.91$\pm$1.42 & 88.3$\pm$48.9            \\
\DM         & 4.56$\pm$0.98 & 254.2$\pm$33.4    & 1.72$\pm$1.19 & 79.1$\pm$42.4            \\
\RL         & 4.91$\pm$0.81 & 197.0$\pm$28.5     & 0.66$\pm$0.72 & 116.2$\pm$27.1            \\
{\bf\CL}   & {\bf0.09$\pm$0.10} & {\bf116.1$\pm$18.0} & {\bf0.58$\pm$0.45} & {\bf39.2$\pm$23.2}\\ \hline
\end{tabular}
\end{table}

\iffalse
\begin{table}[bp]
\scriptsize
\centering
\caption{Generalization results on large instances: Primal gap (PG) (in percent) and primal integral (PI) at 60 minutes runtime cutoff, averaged over 100 instances. ``$\downarrow$'' means the lower the better. \label{res::bigtable60}}

\begin{tabular}{c|rr|rr}
\hline    &  \multicolumn{1}{c}{PG (\%) $\downarrow$}         & \multicolumn{1}{c|}{PI $\downarrow$} & \multicolumn{1}{c}{PG (\%) $\downarrow$} &\multicolumn{1}{c}{PI $\downarrow$} \\ \hline

                       & \multicolumn{2}{c|}{MVC-L}                                                                        & \multicolumn{2}{c}{MIS-L}                                                                         \\ \hline
BnB                    & {2.41$\pm$0.40} & 130.2$\pm$11.1            & {6.29$\pm$1.62} & 285.1$\pm$18.2           \\
\RANDOM                 & {0.38$\pm$0.24} & 22.7$\pm$8.0              & {\bf0.11$\pm$0.08} & 19.0$\pm$3.1             \\
\LBRELAX                 & {0.46$\pm$0.23} & 48.4$\pm$7.5             & {0.91$\pm$0.16} & 68.6$\pm$5.5             \\
\DM               & {0.27$\pm$0.23} & {21.2$\pm$8.1}               & {0.29$\pm$0.15} & { 27.1$\pm$5.5}             \\
\RL               & {0.59$\pm$0.30} & 37.3$\pm$9.6             & {0.14$\pm$0.12} & 18.9$\pm$4.1            \\
{\bf\CL}                & {\bf 0.05$\pm$0.04} & {\bf 9.1$\pm$3.4}             & { 0.12$\pm$0.11} & {\bf12.9$\pm$4.4}            \\ \hline
\multicolumn{1}{l|}{} & \multicolumn{2}{c|}{CA-L}                                                                         & \multicolumn{2}{c}{SC-L}                                                                          \\ \hline
BnB        & 2.74$\pm$1.87 & 320.9$\pm$83.1           & 1.54$\pm$1.33 & 115.0$\pm$42.5            \\
\RANDOM     & 5.37$\pm$0.75 & 229.2$\pm$24.4          & 3.31$\pm$1.79 & 166.4$\pm$61.3         \\
\LBRELAX    & 1.61$\pm$1.50 & 153.0$\pm$50.3          & 1.91$\pm$1.42 & 88.3$\pm$48.9            \\
\DM         & 4.56$\pm$0.98 & 254.2$\pm$33.4    & 1.72$\pm$1.19 & 79.1$\pm$42.4            \\
\RL         & 4.91$\pm$0.81 & 197.0$\pm$28.5     & 0.66$\pm$0.72 & 116.2$\pm$27.1            \\
{\bf\CL}   & {\bf0.09$\pm$0.10} & {\bf116.1$\pm$18.0} & {\bf0.58$\pm$0.45} & {\bf39.2$\pm$23.2}\\ \hline
\end{tabular}
\end{table}

\fi
\begin{figure*}[bpth]
    \centering
    \includegraphics[width=0.8\textwidth]{figure_all/legend_ML_horizontal_timeVSobj_3601.eps}

    \subfloat[MVC-S (left) and MVC-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_survival_rate_0.0100MVC_BA70_1000_3601_nh_50.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_survival_rate_0.0100MVC_BA70_2000_3601_nh_50.eps}
    }
    \subfloat[MIS-S (left) and MIS-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_survival_rate_0.0100INDSET_ER4_6000_3601_nh_3000.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_survival_rate_0.0100INDSET_ER4_12000_3601_nh_6000.eps}
    }\\
    \subfloat[CA-S (left) and CA-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_survival_rate_0.0100CAT_2000_4000_0.7_3601_nh_1000.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_survival_rate_0.0100CAT_4000_8000_0.7_3601_nh_1000.eps}
    }
    \subfloat[SC-S (left) and SC-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_survival_rate_0.0100SC_4000_5000_3601_nh_150.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_survival_rate_0.0100SC_8000_5000_3601_nh_150.eps}
    }
    \caption{The survival rate (the higher the better) over 100 test instances as a function of runtime to meet primal gap threshold 1.00\%. For ML approaches, the policies are trained on only small training instances but tested on both small and large test instances.  \label{res::survival}}
\end{figure*}

\begin{figure*}[btp]
    \centering
    \includegraphics[width=0.8\textwidth]{figure_all/legend_ML_horizontal_timeVSobj_3601.eps}

    \subfloat[MVC-S]
    {
        \includegraphics[width=0.23\textwidth]{figure_all/ML_winRate_MVC_BA70_1000_3601_nh_50.eps}

    }
    \subfloat[MIS-S]
    {
        \includegraphics[width=0.23\textwidth]{figure_all/ML_winRate_INDSET_ER4_6000_3601_nh_3000.eps}

    }
    \subfloat[CA-S]
    {
        \includegraphics[width=0.23\textwidth]{figure_all/ML_winRate_CAT_2000_4000_0.7_3601_nh_1000.eps}

    }
    \subfloat[SC-S]
    {
        \includegraphics[width=0.23\textwidth]{figure_all/ML_winRate_SC_4000_5000_3601_nh_150.eps}

    }
    \caption{The best performing rate (the higher the better) as a function of runtime on 100 small instances (see Appendix for results on large instances).  The sum of the
best performing rates at a given runtime might sum up greater than 1 since ties are counted multiple times.\label{res::winRate}}
\end{figure*}

\iffalse

\begin{figure*}[btp]
    \centering
    \includegraphics[width=0.8\textwidth]{figure_all/legend_ML_horizontal_timeVSobj_3601.eps}

    \subfloat[MVC-S (left) and MVC-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_winRate_MVC_BA70_1000_3601_nh_50.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_winRate_MVC_BA70_2000_3601_nh_50.eps}
    }
    \subfloat[MIS-S (left) and MIS-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_winRate_INDSET_ER4_6000_3601_nh_3000.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_winRate_INDSET_ER4_12000_3601_nh_6000.eps}
    }\\
    \subfloat[CA-S (left) and CA-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_winRate_CAT_2000_4000_0.7_3601_nh_1000.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_winRate_CAT_4000_8000_0.7_3601_nh_1000.eps}
    }
    \subfloat[SC-S (left) and SC-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_winRate_SC_4000_5000_3601_nh_150.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_winRate_SC_8000_5000_3601_nh_150.eps}
    }
    \caption{The best performing rate as a function of time. \label{res::winRate}}
\end{figure*}

\begin{figure*}[btp]
    \centering
    \includegraphics[width=0.8\textwidth]{figure_all/legend_ML_horizontal_timeVSobj_3601.eps}

    \subfloat[MVC-S (left) and MVC-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_virtualBest_MVC_BA70_1000_3601_nh_50.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_virtualBest_MVC_BA70_2000_3601_nh_50.eps}
    }
    \subfloat[MIS-S (left) and MIS-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_virtualBest_INDSET_ER4_6000_3601_nh_3000.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_virtualBest_INDSET_ER4_12000_3601_nh_6000.eps}
    }\\
    \subfloat[CA-S (left) and CA-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_virtualBest_CAT_2000_4000_0.7_3601_nh_1000.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_virtualBest_CAT_4000_8000_0.7_3601_nh_1000.eps}
    }
    \subfloat[SC-S (left) and SC-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_virtualBest_SC_4000_5000_3601_nh_150.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_virtualBest_SC_8000_5000_3601_nh_150.eps}
    }
    \caption{The gap to virtual best as a function of time. \label{res::virtualBest}}
\end{figure*}

\begin{figure*}[tbp]
    \centering
    \includegraphics[width=0.8\textwidth]{figure_all/legend_ML_horizontal_timeVSobj_3601.eps}

    \subfloat[MVC-S (left) and MVC-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSobj_MVC_BA70_1000_3601_nh_50.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSobj_MVC_BA70_2000_3601_nh_50.eps}
    }
    \subfloat[MIS-S (left) and MIS-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSobj_INDSET_ER4_6000_3601_nh_3000.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSobj_INDSET_ER4_12000_3601_nh_6000.eps}
    }\\
    \subfloat[CA-S (left) and CA-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSobj_CAT_2000_4000_0.7_3601_nh_1000.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSobj_CAT_4000_8000_0.7_3601_nh_1000.eps}
    }
    \subfloat[SC-S (left) and SC-L (right).]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSobj_SC_4000_5000_3601_nh_150.eps}
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSobj_SC_8000_5000_3601_nh_150.eps}
    }
    \caption{The primal bound as a function of time, averaged over 100 instances.\label{res_all::bound}}
\end{figure*}
\fi

\paragraph{Hyperparameters}
We conduct experiments on 2.5GHz Intel Xeon Platinum 8259CL CPUs with 32 GB memory. Trainings are done on a NVIDIA A100 GPU with 40 GB memory. All experiments use the hyperparameters described below unless stated otherwise.
We use SCIP (v8.0.1) \cite{BestuzhevaEtal2021OO} to solve the sub-ILP in every iteration of LNS.
To run LNS, we find an initial solution by running SCIP for 10 seconds. We set the time limit to 60 minutes to solve each instance and 2 minutes for solving the sub-ILP in every LNS iteration.
All approaches require a neighborhood size $k^t$ in LNS, except for BnB and \RL ($k^t$ in \RL is defined implicitly by how the policy is used). For \LBRELAX, \DM and \CL, the initial neighborhood size $k^0$ is set to  $100, 3000, 1000$ and $150$ for MVC, MIS, CA and SC, respectively, except $k^0$ is set to $150$ for SC for \DM; for \RANDOM, it is set to $200, 3000, 1500$ and $200$ for MVC, MIS, CA and SC, respectively. All approaches use adaptive neighborhood sizes with  $\gamma=1.02$ and $\beta=0.5$, except for BnB and \RL. For \DM, when applying its learned policies, we use the sampling methods on MVC and CA instances and the greedy method on SC and MIS instances. For \CL, the greedy method is used on all instances.
Additional details on hyperparameter tunings are provided in Appendix.

For data collection, we use different neighborhood sizes $k^0=50, 500, 200$ and $50$ for MVC, MIS, CA and SC, respectively, which we justify in Section \ref{sec::expresults}. We set $\gamma=1$ and run LNS with LB until no new incumbent solution found. The runtime limit for solving LB in every iteration is set to 1 hour.
For training, we use the Adam optimizer \cite{kingma2015adam} with learning rate $10^{-3}$. We use a batch size of 32 and train for 30 epochs (the training typically converges in less than 20 epochs and 24 hours).

\subsection{Results} \label{sec::expresults}

Figure \ref{res::gap} shows the primal gap as a function of runtime.  Table \ref{res::smalltable60}  presents the average primal gap and primal integral at 60 minutes runtime cutoff on small and large instances, respectively (see results at 15, 30 and 45 minutes runtime cutoff in Appendix). Note that we were not able to reproduce the results on CA-S and CA-L reported in \citet{wu2021learning} for \RL despite using their code and repeating training with five random seeds. \CL shows significantly better  anytime performance than all baselines on all problems, achieving the smallest average primal gap and primal integral. It also demonstrates strong generalization performance on large instances unseen during training. Figure \ref{res::survival} shows the survival rate to meet the $1.00\%$ primal gap threshold. \CL achieves the best survival rate at 60 minutes runtime cutoff on all instances, except that, on SC-L, its final survival rate is slightly worse than \RL  but it achieves the rate with much shorter runtime. On MVC-L, MIS-S and MIS-L instances, several baselines achieve the same survival rate as \CL but it always achieves the rates with the shortest runtime. Figure \ref{res::winRate} shows the best performing rate on the small test instances where \CL consistently performs best on  50\% to 100\% of the instances.  In Appendix, we present strong results in comparison with two more baselines and on one more performance metric.

\begin{figure*}[btp]
    \centering
\scriptsize

\begin{tabular}{c|rr|rr|rr|rr}
\hline

                       & \multicolumn{2}{c|}{MVC-S}  & \multicolumn{2}{c|}{MIS-S}  & \multicolumn{2}{c|}{CA-S}  & \multicolumn{2}{c}{SC-S}                                                                         \\ \hline
                                              & \multicolumn{1}{c}{NH size}             & \multicolumn{1}{c|}{Runtime}         & \multicolumn{1}{c}{NH size}             & \multicolumn{1}{c|}{Runtime}& \multicolumn{1}{c}{NH size}             & \multicolumn{1}{c|}{Runtime}& \multicolumn{1}{c}{NH size}             & \multicolumn{1}{c}{Runtime} \\ \hline

LB  &100 & 3600$\pm$0  &  3,000 & 3600$\pm$0 & 1,000 & 3600$\pm$0& 100 &3600$\pm$0\\
LB (data collection) & 50 & 3600$\pm$0 & 500 &3600$\pm$0 &   200 & 3600$\pm$0& 50 &3600$\pm$0\\
\DM & 100 & 2.1$\pm$0.1& 3,000 & 1.3$\pm$0.2 & 1,000 & 20.8$\pm$13.1 & 150 &120.9$\pm$1.3 \\
\CL & 100 & 2.2$\pm$0.1& 3,000 & 1.3$\pm$0.1 & 1,000 & 25.1$\pm$15.3 & 100 & 50.1$\pm$10.4\\ \hline
\end{tabular}
    \includegraphics[width=0.50\textwidth]{figure_all/legend_AS_horizontal_timeVSobj_3597.eps}

    \subfloat[MVC-S]
    {
        \includegraphics[height=2.8cm]{figure_all/AS_primalBoundIteration_MVC_BA70_1000_3597_nh_50_iter11.eps}
    }
    \subfloat[MIS-S]
    {
        \includegraphics[height=2.8cm]{figure_all/AS_primalBoundIteration_INDSET_ER4_6000_3597_nh_3000_iter11.eps}
    }
    \subfloat[CA-S]
    {
        \includegraphics[height=2.8cm]{figure_all/AS_primalBoundIteration_CAT_2000_4000_0.7_3597_nh_1000_iter11.eps}
    }
    \subfloat[SC-S]
    {
        \includegraphics[height=2.8cm]{figure_all/AS_primalBoundIteration_SC_4000_5000_3598_nh_150_iter11.eps}
    }
    \caption{The primal bound (the lower the better) as a function of number of iterations, averaged over 100 small test instances. LB and LB (data collection) are LNS with LB using the neighborhood sizes fune-tunded for \CL and for data collection, respectively. The table shows the neighborhood size (NH size) and the average runtime in seconds  (with standard deviations) per iteration for each approach.   \label{res::boundperIteration}}
\end{figure*}

\paragraph{Comparison with LB (the Expert)}

Both \DM and \CL learn to imitate LB. On the small test instances, we run LB with two different neighborhood sizes, one that is fine-tuned in data collection and the other the same as \CL, for 10 iterations and compare its per iteration performance with \DM and \CL.  This allows us to compare the quality of the learned policies to the expert independently of their speed. The runtime limit per iteration for LB is set to 1 hour. Figure \ref{res::boundperIteration} shows the primal bound as a function of number of iterations. The table in the figure summarizes the neighborhood sizes and the average runtime per iteration. For LB, the result shows that the neighborhood size affects the overall performance. Intuitively, using a larger neighborhood size in LB allows LNS to find better incumbent solutions due to being able to explore larger neighborhoods. However, in practice, LB becomes less efficient in finding good incumbent solutions as the neighborhood size increases, sometimes even performs worse than using a smaller neighborhood size (the one for data collection).  The neighborhood size for data collection is fine-tuned on validation sets to achieve the best primal bound upon convergences, allowing the ML models to observe demonstrations that lead to as good primal bounds as possible in training. However, when using the ML models in testing, we have the incentive to use a larger neighborhood size and fine tune it since we no longer suffer from the bottleneck of LB.  We therefore fine tune the neighborhood sizes for \DM and \CL separately on validation sets. \CL has a strong per-iteration performance that is consistently better than \DM. With the fine-tuned neighborhood size, \CL even outperforms the expert that it learns from (LB for data collection) on MIS-S and CA-S.

\paragraph{Ablation Study}

We evaluate how contrastive learning  and  two enhancements  contribute to \CL's performance. Compared to \DM,  \CL uses (1) addition features from \citet{khalil2016learning} and (2) GAT instead of GCN. We denote by ``FF'' the full feature set  used in \CL and ``PF'' the partial feature set in \DM. In addition to \DM and \CL, we evaluate the performance of \DM with FF and GAT (denoted by IL-LNS-GAT-FF), \CL with GCN and PF (denoted by CL-LNS-GCN-PF) as well as \CL with GAT and PF (denoted by CL-LNS-GAT-PF) on MVC-S and CA-S. Figure \ref{res::ablation_gap} shows the primal gap as a function of runtime. Table \ref{res::ablation_table} presents the primal gap and primal integral at 60 minutes runtime cutoff. The result shows that IL-LNS-GAT-FF, imitation learning with the two enhancements, still performs worse than CL-LNS-GCN-PF without any enhancements.  CL-LNS-GCN-PF and CL-LNS-GAT-PF perform similarly in terms of the primal gaps but CL-LNS-GAT-PF has better primal integrals, showing the benefit of replacing GCN with GAT. On MVC-S, three variants of \CL have similar average primal gaps and on CA-S, \CL has better average primal gap than the other two variants. But adding the two enhancement helps improve the primal integral, leading to the overall best performance of \CL on both MVC-S and CA-S.
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.5\textwidth]{figure_all/legend_AS_2row_timeVSobj_3599.eps}

    \subfloat[MVC-S]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSgap_MVC_BA70_1000_3599_nh_50.eps}
    }
    \subfloat[CA-S]
    {
        \includegraphics[width=0.24\textwidth]{figure_all/ML_timeVSgap_CAT_2000_4000_0.7_3599_nh_1000.eps}
    }

    \caption{Ablation study: The primal gap (the lower the better) as a function of time, averaged over 100 small test instances.\label{res::ablation_gap}}
\end{figure}

\begin{table}[htbp]

\centering
\caption{Ablation study: Primal gap (PG) (in percent) and primal integral (PI) at 60 minutes runtime cutoff, averaged over 100 small test instances and their standard deviations.
``$\downarrow$'' means the lower the better. \label{res::ablation_table}}

\scriptsize
\begin{tabular}{c|rr|rr}
\hline
                       & \multicolumn{1}{c}{PG (\%) $\downarrow$}         & \multicolumn{1}{c|}{PI $\downarrow$}     & \multicolumn{1}{c}{PG (\%) $\downarrow$}         & \multicolumn{1}{c}{PI $\downarrow$} \\ \hline

                       & \multicolumn{2}{c|}{MVC-S}                                                                        & \multicolumn{2}{c}{CA-S}                                                                         \\ \hline
IL-LNS(-GCN-PF)        & {0.29$\pm$0.23} & {19.2$\pm$10.2} &  1.09$\pm$0.51 & 90.0$\pm$20.8            \\
IL-LNS-GAT-FF       & 0.24$\pm$0.17 & 15.3$\pm$7.3 & 1.13$\pm$0.63 & 78.9$\pm$22.7 \\ \hline
CL-LNS-GCN-PF       & 0.17$\pm$0.10 & 11.4$\pm$8,8 & 0.75$\pm$0.40& 57.9$\pm$21.2\\
CL-LNS-GAT-PF       & {\bf0.16$\pm$0.09} & 10.1$\pm$0.6 & 0.76$\pm$0.39 & 53.8$\pm$22.1 \\
{\bf CL-LNS(-GAT-FF)}      & { 0.17$\pm$0.09} & {\bf 8.7$\pm$6.7} &  {\bf 0.65$\pm$0.32} & {\bf50.7$\pm$22.7} \\ \hline

\end{tabular}
\end{table}

\section{Conclusion}

We proposed \CL, that uses a contrastive loss to learn efficient and effective destroy heuristics in LNS for ILPs. We presented a novel data collection process tailored for \CL and used GAT with a richer set of features to further improve its performance. Empirically, \CL significantly outperformed state-of-the-art approaches on four ILP benchmarks w.r.t. to the primal gap, the primal integral, the best performing rate and the survival rate. \CL achieved good generalization performance on out-of-distribution instances.
It is future work to learn policies that can generalize across problem domains. \CL does not guarantee optimality and it is also interesting future work to integrate it in BnB for which many other learning techniques are developed. Our approach is closely related to and could be useful for many problems of identifying substructures in combinatorial searches, for example, identifying backdoor variables in ILPs \cite{ferber2022learning} and selecting neighborhoods in LNS for other COPs.


\section{Introduction} \label{sec:intro}

Since its first indication by Zwicky \citep{Zwicky1933,Zwickly1937}, dark matter has long been thought to be an essential ingredient to explain the cosmic structure formation across a wide range of observations. One important consequence, supported by various observations, is that dark matter forms,  at the macroscopic level, a smooth distribution with virtually null initial local velocity dispersion, referred to as cold dark matter (CDM) \citep{Peebles1982,Peebles1984,Blumenthal_etal1984}.
According to a widely accepted scenario, a collapse of CDM occurs within
a cosmic web, leading to the formation of self-gravitating bound objects called dark matter halos. The late-time evolution of such halos is driven by the continuous accretion of surrounding dark matter onto the halo center and successive mergers with other halos, resulting in the structure of CDM halos exhibiting several characteristic features.
One prominent feature, found in cosmological $N$-body simulations, is the so-called Navarro-Frenk-White (NFW)   profile, first suggested by \citet{Navarro_Frenk_White1996,Navarro_Frenk_White1997}, which has had a significant impact on observations as a testing ground for the CDM paradigm. Another striking feature is the power-law nature of the pseudo-phase-space density profile defined by $Q(r)\equiv\rho(r)/\sigma^3(r)$, with $\rho(r)$ and $\sigma(r)$  being  respectively the radial density and velocity  dispersion profile  \citep{Taylor_Navarro2001,Navarro_etal2010,Ludlow_etal2010}. The slope found in simulations closely matches the  prediction of the Bertschinger’s secondary infall model \citep{Bertschinger1985},
indicating that the structure of halos is built up in a self-similar manner. These properties imply that there is something more fundamental underlying them, which could aid in understanding the origin of universalities.

In this {\it Letter}, we present an alternative way to characterize halo density profiles based on phase-space information. It is worth noting that the formation of CDM halos is accompanied by
a shell crossing at an early phase of matter accretion,
followed by a multi-stream motion of matter distribution. The multi-stream nature of halos has garnered recent attention, highlighted with a renewed interest,
as the outer boundary of the multi-stream region
serves as a natural boundary of halos, and
is clearly demarcated by a radial caustic, manifested as a local steepening of the density profile, referred to as the splashback radius \cite[e.g.,][]{Diemer_Kravtsov2014,Adhikari_Dalal_Chamberlain2014,More_Diemer_Kravtsov2015}.
Motivated by these findings, \cite{Sugiura20} developed a method using an extension of the SPARTA algorithm in \cite{Diemer2017} to reveal the multi-stream nature of halos at the outer regions and they found that about $30$\% of halos are well-described by the self-similar solution of \citet{Fillmore_Gorldreich1984}. In this {\it Letter}, by substantially refining their analysis based on high-resolution simulations with finely sampled snapshots out to an early
halo formation, we are able to unveil the innermost parts of the multi-stream region, where we find that halos exhibit a universal feature in each multi-stream distribution.

\section{Method} \label{sec:method}

We analyze cosmological $N$-body simulations performed in a flat $\Lambda$CDM cosmology, which is consistent with recent observations of cosmic microwave background radiation~\citep{Planck2015cosmology}. We mainly analyze the simulation that follows the movements of $500^3$ particles in a comoving box with a side length of $41\,h^{-1}\mathrm{Mpc}$ using the Tree Particle-Mesh code \textsc{Ginkaku}~(Nishimichi, Tanaka \& Yoshikawa, in preperation). We employ a softening length of $4.1\,h^{-1}\mathrm{kpc}$, which we denote by $r_\mathrm{LR}$ in what follows. The snapshots of the particles are saved at $1,001$ redshifts, evenly spaced between $z=0$ and $5$, providing a dense sampling to accurately determine the number of apocenter passages (denoted by $p$ in what follows) up to $\sim 50$, following the method of \citet{Sugiura20} with minor modifications.

We first select relaxed halos from those identified by \textsc{Rockstar}~\citep{Rockstar} at $z=0$, by imposing a cut in the spin parameter and the offset between the center of mass and the density peak~\citep{Klypin16}. We also discard subhalos according to the consistency between the exact spherical-overdensity mass and that listed in the Rockstar catalog. We then trace the main progenitor by following the particles within the virial radius back in time, updating the center and the list of member particles using the shrinking-sphere method at each snapshot until we reach the first snapshot at $z=5$ or the number of member particles falls below $1,000$. Our final halo trajectories are defined as the center of mass of the $1,000$ fixed member particles, which are closest in phase space to the center of the main progenitor at the highest redshift to which we can trace the progenitor with at least $1,000$ particles. We next follow forward in time the center of mass of these fixed particles to obtain a smooth trajectory robust to merger events. We
monitor the velocities and positions of all surrounding particles that are within $2.5 \, R_\mathrm{vir}$ at $z=0$ relative to the center of the progenitor. We
define and count the apocenter passage
for each particle when the relative velocity changes from outgoing to infalling and the relative position has orbited at least $90^\circ$ from the previous apocenter passage~\citep{Sugiura20}. These specific choices are found to be robust for the determination of the number of apocenter passages up to $\sim 50$.

\begin{figure}[t!]
\epsscale{1.2}
\plotone{graphs/radial_properties2_halorank1_snapnum1000.pdf}
\caption{Radial density profile (upper) and phase-space distribution (lower) of a halo with $M_\mathrm{vir}=1.49\times 10^{14}\,h^{-1}M_\odot$. The upper panel shows the decomposition of the total density profile (highest line) into the contributions from $N$-body particles with different numbers of apocenter passages, represented by colors ranging from $p=1$ (dark blue) to $p=50$ (dark red) . The lower panel displays the distribution of individual particles, with the same color coding. The infalling component, $p=0$, is depicted in gray.
\label{fig:indivhalo}}
\end{figure}

In Fig.~\ref{fig:indivhalo}, we present the radial density profile and phase-space distribution of a representative halo with mass $M_{\rm vir}=1.49\times10^{14}\,h^{-1}\,M_\odot$, color-coded by the number of apocenter passages, $p$. It is apparent that particles with a high value of $p$ tend to be concentrated at smaller radii, leading to an increase in density
and a reduction in velocity dispersion, resulting in an onion-like multi-stream structure in the phase-space distribution. On the other hand, the density profiles exhibit similar features, with the inner and outer slopes converging to a specific value regardless of $p$. In the following sections, we will further analyze this behavior for halos with different properties.

In order to study the convergence, we have conducted a higher-resolution simulation with $2,000^3$ particles with an identical initial Gaussian random field. However, storing as many as $\gtrsim 1,000$ snapshots from this simulation requires a significant amount of disk space, and
an accurate apocenter count would be costly. Therefore, we only use this run to verify the density profile at $z=0$.
In the following discussion, we refer to this simulation as HR, while the one with $500^3$ particles is called LR.
The softening scale for HR is $r_\mathrm{HR} = 1.025\,h^{-1}\mathrm{kpc}$.

\section{Results} \label{sec:results}

In order to systematically and quantitatively investigate the radial density profile of particles with $p$ apocenter passages, we divide the halos into four mass bins (Table \ref{tab:halo_mass_bin}). For each mass bin, we measure the stacked profile for each stream by rescaling the radial coordinate by the virial radius and the density by the virial mass for individual halos.
We
find that the stacked profile for each $p$, as presented in the upper panels and the lower left panel in Figure~\ref{fig:general}, are well-fitted by the following functional form\footnote{
The fitting analysis employs the standard $\chi^2$ method, with weights determined by the inverse variance among the stacked halos. Radial bins with $r > 1.2\, r_\mathrm{LR}$ are considered.}
:
\begin{equation}
        \rho_{\rm stream}(r;p)=\frac{A(p)}{\{r/S(p)\}
        \bigl[1+\{r/S(p)\}^7\bigr]},
    \label{eq:double_powerlaw_fit}
\end{equation}
where the characteristic scale $S(p)$ and density $A(p)$ are given as a function of $p$. Due to the rescaling in the stacking, the functions, $S(p)$ and $A(p)$, are quite similar among the four mass bins, as shown in Figure~\ref{fig:fitparams}.
We find that these functions are well captured by
\begin{equation}\label{eq:fitA}
\begin{split}
    \log_{10}{\Bigl\{A_\mathrm{fit}(p)/\overline{\rho}_{\rm m}\Bigr\}}  = 4.89-0.119\log_{10}{\left(M_\mathrm{vir,10}\right)}  \\ +\Bigl\{-3.89+0.243\log_{10}{\left(M_\mathrm{vir,10}\right)}\Bigr\}\,p^{-9/40},
\end{split}
\end{equation}
\begin{align}\label{eq:fitS}
\begin{split}
    \log_{10}{\Bigl\{S_\mathrm{fit}(p)/R_{\rm vir}\Bigr\}} =
    2.46-0.0474\log_{10}{\left(M_\mathrm{vir,10}\right)}  \\ +
    \Bigl\{-2.29-0.0639\log_{10}{\left(M_\mathrm{vir,10}\right)}\Bigr\}\,p^{1/8},
\end{split}
\end{align}
including the weak mass dependence,
where $M_\mathrm{vir,10}$ is defined by $M_\mathrm{vir}/10^{10}h^{-1}M_\odot$.

In Figure~\ref{fig:general}, the stacked profiles are presented up to $p=40$.
The horizontal axis is scaled by $(p/R_\mathrm{vir})$ and the vertical axis by $(10^{p/4}\,r/\overline{\rho}_{\rm m})$ for improved visibility.
The errorbars indicate the standard error of the stacked profiles.
The profiles are in agreement with Equation~(\ref{eq:double_powerlaw_fit}) over a wide range of radii and mass scales. A sharp cutoff is observed in the profile at large radii, which is consistent with our model with the asymptotic slope of $-8$.
This is in line with the findings of \citet{2022arXiv220503420D}, who identified an exponential form with radial-dependent exponents for the profiles of ``orbiting'' particles, after the first pericenter passage.
More notably, the inner slope tends to be consistent with $-1$ for most cases, except for orbits with $p\lesssim10$ for the XL sample, which exhibit a shallower slope.
This is likely due to the sensitivity of these low-$p$ orbits to recent mass accretion or merger history~\citep[e.g., ][]{Sugiura20}. However, this trend tends to be erased after several orbits, reaching a universal slope for $p\gtrsim10$, indicating a self-similar growth of phase-space structure.

To quantitatively assess the double power-law nature of each stream, we
compare the total density profile from HR for
halos that have been matched with LR to the prediction obtained by summing the individual double power-law profiles described by Equations~(\ref{eq:double_powerlaw_fit})--(\ref{eq:fitS})\footnote{In the plot, the summation is conservatively taken up to $p=3,000$. The change in density is less than $0.2\%$ over the plotted range when we instead stop at $p=300$.}.
The results are shown in Figure \ref{fig:fitsum}, where the solid lines with shaded regions represent the prediction based on the double power-law model, taking into account the uncertainties in the numerical coefficients in the fit.
Our model is in good agreement with HR for all four mass bins. Notably, we can recover the profile below $r_\mathrm{LR}$, despite the fact that the individual profiles for each $p$ are fitted above $1.2\,r_\mathrm{LR}$ and only up to $p=40$.
This suggests that the model effectively extrapolates the mass distribution to large values of $p$ beyond the resolution limit of LR. In the lower panel, we can observe the transition of the slope from $-3$ to $-1$ in different models\footnote{The logarithmic slope is estimated from discrete simulation data points with statistical noise using \texttt{GEORGE} python package~\citep{GEORGE} for Gaussian Process.}.
\begin{figure*}[ht!]
\includegraphics[width=\linewidth]{graphs/rdpfit_rhor_all_APJ_snapnum1000.pdf}
\caption{Stacked radial density profiles of $N$-body particles with even number of apocenter passages, ranging from $p=4$ to $40$.
The four mass bins are displayed in the upper three and the lower left panel for S, M, L and XL, respectively.
Additionally, the lower middle and right panels show the results obtained from $460$ halos in the mass range $[4.10\time10^{11},\,2.30\times10^{12}]\,h^{-1}M_\odot$,
which are further divided
into the two subsamples based on the concentration parameter $c_{\rm vir}$ and accretion rate $\Gamma_{\rm dyn}$, respectively (see text in detail). In each panel, the fitted results with Equation~(\ref{eq:double_powerlaw_fit}) are depicted as solid lines.}
\label{fig:general}
\end{figure*}

\begin{figure}[ht!]
\epsscale{1.15}
\plotone{graphs/stacked_dpfitparams_A10B80_snapnum1000.pdf}
\caption{
Dependence of the characteristic density $A$ (upper) and scale $S$ (lower) on the number of apocenter passages, $p$, as determined by fitting to Equation~(\ref{eq:double_powerlaw_fit}) in different symbols four mass bins (see legend). The thin solid curves represent the fitting formulae, Eqs.~(\ref{eq:fitA}) and (\ref{eq:fitS}). For comparison, predictions of the Fillmore-Goldreich self-similar solutions are also shown,
for specific values of the parameter
$\epsilon$ ($1/15$, $1/6$ and $1$).
In plotting these predictions, we identify the position of radial caustics in the self-similar solutions with the characteristic scale $S(p)$, and derive $A(p)$ by equating the masses contained in each stream. The shaded regions for the predictions indicate uncertainty in identifying $S(p)$ with the position of the $p$-th or $(p+1)$-th radial caustics of the self-similar solutions.}
\label{fig:fitparams}
\end{figure}

\begin{figure*}[ht!]
\epsscale{1.15}
\plotone{graphs/stacked_dpfitsum_A10B80_snapnum1000_2.pdf}
\caption{
Comparison of the total density profile between our model
($\sum\rho_{\rm stream}$) and HR ($\rho_\mathrm{HR}$). The upper and lower panels respectively show the fractional difference with respect to HR, i.e., $(\rho-\rho_{\rm HR})/\rho_{\rm HR}$, and the logarithmic slope, $d\log\rho/d\log r$.
The results for four mass bins are presented separately in each panel for scales above twice the softening length of HR, $r_\mathrm{HR}$.
The shaded regions indicate the estimated uncertainties in the prediction, which are propagated from the statistical error in the stacked profile through the uncertainties in the fitting parameters.
We also plot the NFW (dashed) and Einasto (dot-dashed; \protect\citet{einasto1965}) profiles, obtained by fitting HR in the range $2 \, r_{\rm HR}\leq r\leq 0.9 \, R_{\rm vir}$.
In the upper panels, the results obtained from a partial summation of the double power-law profile up to $p=40$ are also plotted (dotted).
The vertical arrows indicate $r_{\rm LR}$, the softening length for LR.}
\label{fig:fitsum}
\end{figure*}

\begin{table}[t!]
  \caption{Halo samples in the LR simulation. The second, third, and fourth columns respectively show the range of halo masses, virial radii, and the number of halos.}
  \label{tab:halo_mass_bin}
  \centering
  \begin{tabular}{cccc}
  \hline\hline
    Sample &
    $10^{11}\,M_\mathrm{vir}\,[h^{-1}M_\odot]$
    &

    $R_\mathrm{vir}\,[h^{-1}\mathrm{Mpc}]$
    & \# of halos \\
    \hline

    S &  $[3.16,\;5.71]$ & $[0.14,\;0.17]$& 300
    \\

    M & $[5.71,\;24.2]$ &$[0.17,\;0.27]$& 300
    \\

    L & $[24.2,\;134]$ &$[0.27,\;0.48]$& 70
    \\

    XL & $[134,\;1530]$ &$[0.48,\;1.08]$& 13
    \\
    \hline
      \end{tabular}
\end{table}

\section{Discussion} \label{sec:discussion}

To further examine the robustness of our results, we analyze a subset of $460$ halos within a specific mass range $[4.10\times10^{11},\,2.39\times10^{12}]\,h^{-1}M_\odot$.
These halos are divided into two sub-samples based on two different criteria.
We employ the concentration parameter $c_{\rm vir}$, which is obtained by fitting the density profiles with the NFW profile for $2\,r_\mathrm{LR}\leq r \leq 0.9\,R_\mathrm{vir}$,
and the mass accretion rate defined by $\Gamma_{\rm dyn}(t)\equiv \{\log[M(t)]-\log[M(t-t_{\rm dyn})]\}/\{ \log[a(t)]-\log[a(t-t_{\rm dyn})]\}$ with $t_{\rm dyn}$ being the dynamical time estimated from halo masses \citep{Diemer2017}\footnote{
We use the virial mass, $M_\mathrm{vir}$, to measure $\Gamma_\mathrm{dyn}$, whereas \cite{Diemer2017} uses $M_\mathrm{200m}$.}.
We divide the halos into two halves, one with high values of these indicators and the other with low values.

The middle bottom (right bottom) panel of Figure \ref{fig:general} depicts the results for two subsamples having low and high values of $c_{\rm vir}$ ($\Gamma_{\rm dyn}$), represented by red and black colors, respectively. Again, a good agreement between the double power-law function and measured profiles is observed over a wide range of $p$. A close look at each stream profile reveals that halos with high concentration or low accretion-rate tend to have a large amplitude $A(p)$ and a large characteristic scale $S(p)$. These trends are particularly evident for larger $p$, suggesting that the universal double power-law feature is established in a self-regulated manner during the orbital motion in the multi-stream region, where the diversity of mass accretion and merger histories tend to be erased and only be imprinted in $A(p)$ and $S(p)$.

In order to gain a deeper understanding of the results obtained in this {\it Letter} from a dynamical viewpoint, it would be beneficial to compare them with self-similar solutions.
While self-similar solutions are only valid in the Einstein-de Sitter universe, the secondary infall model of \cite{Bertschinger1985} has been shown to reproduce the pseudo phase-space density of $Q(r)\propto r^{-1.875}$ found in simulations in the $\Lambda$CDM model.
One such solution that can be used to describe the multi-stream flow in phase space  is the spherically symmetric solutions put forth by \citet{Fillmore_Gorldreich1984}.
Under a stationary matter accretion, this model includes the Bertschinger's secondary infall model as a special case.
Recent work by \citet{Sugiura20} has made a direct comparison of these predictions with radial multi-stream structures obtained from simulations up to $p=5$. Identifying the position of radial caustics in self-similar solutions with the characteristic scale of the double power-law profile in Equation~(\ref{eq:double_powerlaw_fit}), it is possible to make predictions for both $A(p)$ and $S(p)$.

In Figure \ref{fig:fitparams}, we compare the predictions of self-similar solutions with our $N$-body results for three values of the model parameter $\epsilon$, which describes the power-law slope of initial density contrast.
Note that the parameter $\epsilon$ is restricted to the range $[0,\,1]$, and the solution with $\epsilon=1$ corresponds to Bertschinger's secondary infall model. Figure \ref{fig:fitparams} shows that none of the solutions consistently explain the trends in both $A(p)$ and $S(p)$, although setting the parameter $\epsilon$ to $1/6$ reproduces the characteristic scale $S(p)$ reasonably well.
The main reason for this failure is that for each stream, the Fillmore-Goldreich solutions predict a
steep inner profile with a logarithmic slope of around $-2$ irrespective of the value of $\epsilon$.
One possible explanation for the shallow inner cusps found in simulations is to introduce the non-zero angular momentum, which can reduce the steepness of the profile near the halo center~\citep{2001MNRAS.325.1397N,2010PhRvD..82j4044Z}.
However, existing solutions allow for the introduction of angular momentum in a very specific manner, and without a broad angular momentum distribution, they fail to describe the shallow inner cusp seen in the profile for each $p$.

We thus conclude that a more comprehensive theoretical study is needed to fully understand the universal features found in this \textit{Letter}, taking into account the complexities associated with mass accretion and merger history.
This may involve exploring the angular momentum distribution or relaxing the symmetry assumptions ~\citep[see][for the latter aspect]{1993ApJ...418....4R,Lithwick_2011}.

\section{Conclusion} \label{sec:conclusion}
In this \textit{Letter}, we have investigated the multi-stream radial structures of dark matter halos in cosmological $N$-body simulations. Our focus is on the radial distribution of dark matter particles within the splashback radius. We use the method developed by \cite{Sugiura20} to trace the trajectories of dark matter particles and quantify the density profile for each stream, which we label by $p$. With the help of $1,001$ snapshots between $z=0$ and $5$, we are able to resolve the multi-stream structure in phase space up to $p=40$.
The radial density profiles for each stream are accurately described by a double power-law function (Equation~\ref{eq:double_powerlaw_fit}), with characteristic density $A(p)$ and scale $S(p)$ well-fitted respectively to Equations ~\eqref{eq:fitA} and ~\eqref{eq:fitS}. These results are consistent across different sample selections based on the concentration parameter and mass accretion rate. We can recover the total density profile by summing up the individual contribution modeled by Equation ~(\ref{eq:double_powerlaw_fit}), which provides a prediction comparable to or even better than the Einasto profile. Our findings suggest that the double power-law nature seen in the stream profiles is universal.
To gain a deeper understanding of these results, we compare them with predictions based on self-similar solutions. We find that the Fillmore-Goldreich solutions (nor Bertschinger's solution as a special case) cannot consistently explain both $A(p)$ and $S(p)$. This suggests that a more comprehensive theoretical study is necessary, taking into account the dynamical complexities associated with halo accretion and merging history.

Finally, the universal features of halos found in this \textit{Letter} are a direct consequence of the cold nature of dark matter, and serve as valuable insights into the physical properties of CDM halos.
While this study has utilized $N$-body simulations and investigated the inner multi-stream structure up to $p=40$,
recent developments in simulating collisionless self-gravitating systems through
Vlasov-Poisson equations offer a promising way to further probe the phase-space structure ~\citep{Yoshikawa_Yoshida_Umemura2013, Hahn_Angulo2016, Sousbie_Colombi2016}. This would provide
a deeper understanding of the physics behind the universal features. Another point worth
further investigating is to scrutinize the radial phase-space structures for alternative dark matter models, as the nature of dark matter
has a significant impact on small-scale structure formation~\citep[e.g.,][for a review]{Bullock_Boylan-Kolchin2017}.
Our method to reveal multi-stream structures can be straightforwardly applied to simulations of other dark matter models.
Any difference in the radial multi-stream structures could provide valuable observational
probes to clarify the nature of dark matter.

\begin{acknowledgments}
We thank St\'ephane Colombi and Takashi Hiramatsu for insightful suggestions and discussions, Shogo Ishikawa and Satoshi Tanaka for comments and discussions. This work was
supported in part by MEXT/JSPS KAKENHI Grant
Number JP19H00677 (TN), JP20H05861, JP21H01081
(AT and TN), and JP22K03634 (TN). We also acknowledge financial support from Japan Science and Technology Agency (JST) AIP Acceleration Research Grant Number JP20317829
(AT and TN).
YE is also supported by JST, the establishment of university fellowships towards the creation of science technology innovation, Grant Number JPMJFS2123.
Numerical computations were carried out
at Yukawa Institute Computer Facility, and Cray XC50 at
Center for Computational Astrophysics, National Astronomical Observatory of Japan.
\end{acknowledgments}