Purity Law for Neural Routing Problem Solvers with Enhanced Generalizability

Thank you very much for your time and effort in reviewing our work. We feel very glad to know you acknowledge our work is intuitive and original, find the statistical analysis is thoroughly conducted, and believe that our paper is well written, logically structured, and presents an interesting perspective.

We will address each of your concerns in detail as follow.

--Q1: Similar law for non-Euclidean problems?

Thank you for raising this important concern. In non-Euclidean problems, there is no longer direct connectivity between vertices, and the metric is given a priori by the edge cost matrix. Given $G=(V, E)$, we define purity order $K_p(e_{ij})$ in non-Euclidean routing problems as follow:

$K _ p(e _ {ij}) = \frac{max _{m _ i \in M _ i} |\xi _ {im _ i}| + max _ {m _ j \in M _ j} |\xi _ {jm _ j}|}{2}$ , where $M _ i$ consists of all $m _i = \arg\max _ {m, l(\xi _ {im}) \leq l(e _ {ij})/2} l(\xi _ {im})$,

and $\xi_{im}$ is the shortest path from vertex $v_i$ to $v_m$.

To understand it intuitively, we first need to find the shortest path starting from the two endpoints, which does not exceed half the length of the edge and has the largest length, and then define the purity order of the two end points as the average of the maximum number of nodes in the shortest path. We take an instance of the non-Euclidean TSP dataset of size 20 in [1] as an example and calculate the pure order according to the above definition. The purity order distribution is as follows:

We can observe that the purity order still exhibits a negative exponential law on non-Euclidean TSP instances. In the future, we will calculate the purity order distribution of optimal solutions for more non-Euclidean TSP instances and add them to the appendix of this paper.

--Q2: More theoretical discussion on purity law.

Thank you for this suggestion. Regarding your first point, we also find the stability of the purity law under different distributions interesting. We speculate that certain types of vertex distributions in instances may satisfy persistent homology in the field of topology, under measured by the "purity order on the optimal solution." However, we have not yet achieved a formal theoretical proof. In future work, we plan to explore more theoretical proofs of this relationship from the perspective of topological data analysis.

Your second point is also interesting and inspiring. Taking the TSP as an example, for instances consisting of regular N-gon endpoints, the negative exponential law is not strictly satisfied. Because the purity order of the edges in the optimal solution is all 0. However, this can be viewed as a limiting case of the negative exponential law. Furthermore, we speculate that there should be special instances that theoretically violate the negative exponential law. In future work, we plan to theoretically construct such special instances or generate them through learning-based methods.

Furthermore, we have rigorously proved an upper bound on the probability that an edge of purity order $k$ is located on the optimal solution under certain conditions. If you are interested, please refer to our response to Reviewer 8oG4.

--Q3: Incorporating purity orders as instance input representation.

Thank you for this suggestion. On one hand, purity order can be incorporated as an edge feature into graph representation learning, such as Graph Convolutional Neural Networks architecture. On the other hand, purity order can be generalized from edges to certain vertex attributes to reflect the purity information of nodes in instances. In future work, we will explore more methods for incorporating purity information into representation learning.

--Q4: Discussion on PUPO combined with heavy-decoder paradigm.

Thank you for this suggestion. The LEHD and BQ-NCO you mentioned are both excellent works, which use imitation learning for training. The loss function of imitation learning is as follows: $l = - \frac{1}{N} \sum_i (y_i \mathop{\log} p_i + (1-y_i) \mathop{\log} (1 - p_i) )$, where $y_i$ is optimal label.

Following the PUPO in our paper, we propose the modified loss function as follow: $\hat{l} = - \frac{1}{N} \sum_i (y_i C(U_{i-1}, \tau_i) \mathop{\log} p_i + (1-y_i) C(U_{i-1}, \tau_i) \mathop{\log} (1 - p_i) ).$

We then compare the two training methods using LEHD as an example, testing them on a few instances of TSPLib. The results are shown below.

As can be seen, PUPO-training improves the model's generalization ability on larger scales, but may exhibit a trade-off on smaller instances. It's worth noting that due to limited time, neither method was fully trained. We believe that thorough training would further improve the performance of both models. A performance comparison of the two models after full training will be supplemented in the appendix.

--Q5: More specific details on training.

Thank you for raising this concern. In this paper, all models are trained on a uniform distribution and then directly tested on the U, C, E, and I distributions. The test results in this paper reflect the performance changes caused by the distribution shift you mentioned. From Table 1 and Table 2 in this paper, PUPO-training improves the model's generalization ability to distribution shift and scale increase.

Reference:

[1] Matrix Encoding Networks for Neural Combinatorial Optimization, NeurIPS 2021

Thank you very much for your time and effort in reviewing our work. We feel very glad to know you acknowledge our purity is novel and impressive, find the theoretical analysis is solid, and believe that our paper is well written.

We will address each of your concerns in detail as follow.

--Q1: Including concorde in the experiments.

Thank you for raising this concern. The test datasets used in this paper is from [1]. In the datasets with multi-distribution, TSP-100 is solved by Gurobi, an exact solver, so the solutions are guaranteed to be truly optimal. TSP-1000, TSP-5000 and TSP-10000 are solved by a SOTA heuristic algorithm LKH3, which can output near-optimal solutions. To settle down your concerns, we use concorde to resolve three large datasets. We report the gap between LKH3 and concorde on TSP-1000 as follow, where negative gap means concorde performs better than LKH3.

It can be observed that LKH3 can reach the optima very nearly. However, due to unacceptable time consumption (one instance need more than 20 hours), concorde fails on the TSP-5000 and TSP-10000. In the future, we will make the concorde-solved optimal values and solutions of each of 200 instances with four distributions in TSP-1000 all open source.

--Q2: Deeper explanation of the mechanism of implicit regularization.

Thank you for this concern. The slight performance decline on small-scale problems derives from the underlying mechanism of implicit regularization of PUPO. Shown in Table 5, we can observe that the PUPO training reduces the sum of the parameter Frobenius norms. It results in an expected but reasonable trade-off on small, identically distributed datasets.

When the scale and distribution of training data is limited and completely same with the test data, PUPO-trained model may underfit the specific patterns and noise of the training set. This is because regularization prevents the model from excessively memorizing "details" of the training data, which might be important for decision-making in a small, same-distributed setting. Vanilla-trained models are more likely to achieve near-optimal fit to the training set in this "ideal" small, same-distributed setting.

However, during vanilla training, tour length serves as the sole reward signal. Lacking guidance from generalizable structure, the model often overfits to these specific patterns in the small training set, resulting in weak generalization. PUPO introduces purity information into the model training process, allowing the model to focus on generalizable patterns even in the small training set, significantly improving cross-distribution and larger-scale generalizability, which is potentially more challenging and relevant in the real world.

On multiple large-scale datasets (1000, 5000, and 10000), our algorithm achieves significant and consistent enhancements compared to vanilla training. In particular, as shown in Table 2, the PUPO-trained ELG-50 achieves a 68.58% improvement on the E-5000 of CVRP. On the distribution-shifted test sets (C, E, I, TSPLib and CVRPLib), the PUPO-trained model also achieves significant cross-distribution generalization, demonstrating the effectiveness of the implicit regularization introduced by PUPO in addressing data distribution-shifted generalization.

In future work, to balance the model's performance on small datasets, we will introduce a tuning parameter $\alpha \geq 0$ and generalize the Purity Weightings to $\hat{W}(U_t, \tau_{t+1}) = W^{\alpha}(U_t, \tau_{t+1})$ to adjust the implicit regularization degree. When $\alpha$ reaches 0, PUPO degenerates to a vanilla policy gradient.

--Q3: More discussion on supervision-learning-based neural solvers. Thank you for this suggestion. For the supervision-learning-based neural solvers, their loss function of imitation learning is as follows: $l = - \frac{1}{N} \sum_i (y_i \mathop{\log} p_i + (1-y_i) \mathop{\log} (1 - p_i) )$, where $y_i$ is optimal label.

Following the PUPO in our paper, we propose the modified loss function as follow: $\hat{l} = - \frac{1}{N} \sum_i (y_i C(U_{i-1}, \tau_i) \mathop{\log} p_i + (1-y_i) C(U_{i-1}, \tau_i) \mathop{\log} (1 - p_i) ).$

We then compare the two training methods using LEHD [2] as an example, testing them on a few instances of TSPLib. The results are shown below.

As can be seen, PUPO-training improves the model's generalization ability on larger scales, but may exhibit a trade-off on smaller instances. It's worth noting that due to limited time, neither method was fully trained. We believe that thorough training would further improve the performance of both models. A performance comparison of the two models after full training will be supplemented in the appendix.

Reference:

[1] INViT: A Generalizable Routing Problem Solver with Invariant Nested View Transformer, ICML 2024

[2] Neural Combinatorial Optimization with Heavy Decoder: Toward Large Scale Generalization, NeurIPS 23

Thank you very much for your time and effort in reviewing our work. We feel glad to know you acknowledge our work is novel and effective, and believe our presenting and motivation are clear.

We will address each of your concerns in detail as follow.

--Q1: Inconsistence between experimental data in this paper with those in related works?

Thanks for your thorough review. In our experiments, to ensure fairness, we conduct each training process for the model under the same hardware conditions and parameters (except for the learning rate). In the original INViT paper, all numerical experiments are performed on an NVIDIA RTX 4090. While in the original ELG paper, models are trained on NVIDIA RTX 6000 Ada with 48GB memory, tested on one NVIDIA RTX 4090. Whereas our experiments are implemented on an NVIDIA RTX 3090. Additionally, the selection of the starting point for the tour is random. These factors may account for the slight differences observed when compared to the results in related works. If the paper is accepted, we will make all the code and detailed training logs open source.

--Q2: Performance decline on small-scale problems.

Thanks for your thorough review. The slight performance decline on small-scale problems derives from the underlying mechanism of implicit regularization of PUPO. Shown in Table 5, we can observe that the PUPO training reduces the sum of the parameter Frobenius norms. It results in an expected but reasonable trade-off on small, identically distributed datasets.

When the scale and distribution of training data is limited and completely same with the test data, PUPO-trained model may underfit the specific patterns and noise of the training set. This is because regularization prevents the model from excessively memorizing "details" of the training data, which might be important for decision-making in a small, same-distributed setting. Vanilla-trained models are more likely to achieve near-optimal fit to the training set in this "ideal" small, same-distributed setting.

However, during vanilla training, tour length serves as the sole reward signal. Lacking guidance from generalizable structure, the model often overfits to these specific patterns in the small training set, resulting in weak generalization. PUPO introduces purity information into the model training process, allowing the model to focus on generalizable patterns even in the small training set, significantly improving cross-distribution and larger-scale generalizability, which is potentially more challenging and relevant in the real world.

On multiple large-scale datasets (1000, 5000, and 10000), our algorithm achieves significant and consistent enhancements compared to vanilla training. In particular, as shown in Table 2, the PUPO-trained ELG-50 achieves a 68.58% improvement on the E-5000 of CVRP. On the distribution-shifted test sets (C, E, I, TSPLib and CVRPLib), the PUPO-trained model also achieves significant cross-distribution generalization, demonstrating the effectiveness of the implicit regularization introduced by PUPO in addressing data distribution-shifted generalization.

In future work, to balance the model's performance on small datasets, we will introduce a tuning parameter $\alpha \geq 0$ and generalize the Purity Weightings to $\hat{W}(U_t, \tau_{t+1}) = W^{\alpha}(U_t, \tau_{t+1})$ to adjust the implicit regularization degree. When $\alpha$ reaches 0, PUPO degenerates to a vanilla policy gradient.

--Q3: More discussion related to algorithm design.

Thank you for this suggestion. We supplemented the theoretical analysis on optimization error of PUPO. First, we denote $\nabla V^{p_{\theta}} = E_{p_{\theta}(\tau | \mathcal{X})} \Biggl[ \sum \limits_{t=2}^{N} \left( L(\tau) - b(X) \right) W(U_{t-1}, \tau_{t}) \nabla \log p_{\theta}(\tau_{t} | \tau_{1:t-1}, X) \Biggr]$,

and $\widehat{\nabla V^{p_{\theta}}} = \sum \limits_{t=2}^{N} \left( L(\tau) - b(X) \right) W(U_{t-1}, \tau_{t}) \nabla \log p_{\theta}(\tau_{t} | \tau_{1:t-1}, X)$.

The update rule of gradient ascent algorithm is $\theta_{t+1} = \theta_t + \eta_t \widehat{\nabla V^{p_{\theta}}}$.

Let us say a function $f$ is $\beta$-smooth if $||\nabla f(\theta)-\nabla f(\theta')|| \leq \beta||\theta-\theta'||$, where the norm $||\cdot||$ is the Euclidean norm. Assume that for all $\theta$, $V^{p_{\theta}}$ is $\beta$-smooth and bounded $V_{\star}$. Suppose the variance is bounded as follows: $E [ ||\widehat{\nabla V^{p_{\theta}}}-\nabla V^{p_{\theta}}||^2 ] \leq \sigma^2$. For $t \leq \beta (V^*-V^{(0)})/\sigma^2$, suppose we use a constant stepsize of $\eta_t = 1/\beta$, and thereafter, we use $\eta_t = \sqrt{2/(\beta T)}$. For all $T$, we can obtain:

$\min\limits_{t \leq T} E[||\nabla V^{(t)}||^2] \leq \frac{2 \beta (V^*-V^{(0)})}{T} + \sqrt{\frac{2\sigma^2}{T}}$

The proof method is similar to [1]. And this optimization error theorem confirms that our proposed algorithm can converge to stationary points.

However, due to the black-box nature of neural networks, proving the approximation ratio and generalization error of learning methods is extremely difficult. The gap evaluation metric in this paper can experimentally reflect the function of the approximation ratio, showing the gap ratio between model performance and the optimal solution. In future work, we hope to provide more comprehensive algorithm design analysis from a theoretical perspective. The above optimization error analysis will be added to the appendix of this paper.

--Q4: Guidance role of supermodularity property to algorithm design.

Thank you for raising this concern. In the late stages of the neural constructive solver's decision-making process, the set of unvisited nodes is quite narrow. Because the agent's early stages occupy underlying better future node selection options, it is likely to choose poor ones that are detrimental to generalization. However, in vanilla training, tour length is the only reward signal, which fails to identify this impact of late-stage decisions. These late-stage decisions, made within a narrow feasible set, are crucial for generalization. The increasing marginal returns of the supermodular function, Purity Availability, effectively capture this tail effect.

Inspired by this motivation, we constructed a Purity Availability metric with supermodular properties to measure the purity potential of the unvisited vertex set. On one hand, since Purity Availability represents the average minimum available purity order, higher values indicate that the state will not be able to generate highly pure structures in the future. For example, in the Figure 5 in appendix, the Purity Availability of strategies (a) and (b) at the current state are 1 and 0, respectively. It means that policy (b) possess ability to construct purer structures in the future. On the other hand, constrained by the necessity to select nodes only from the unvisited vertex set, later decision stages often lead to states that are unfavorable for generalization and difficult to distinguish under tour length metrics.

Due to supermodularity, the marginal benefit of Purity Availability, $\phi (U_t) - \phi (U_{t-1})$, is greater at later stages, making the purity measure in the later decision phases more significant. Inspired by supermodularity, we designed the purity cost based on the marginal benefit of state and the purity order of the actions. Furthermore, we designed purity weights based on the concept of future discounting to comprehensively assess the future purity potential of intermediate states during policy transitions.

--Q5: Incorporating purity cost directly into the reward.

Thank you for this suggestion. We use the sum of the pure orders of the edges in the tour as a regularization term and design the following modified loss function: $\hat{L}(\tau) = L(\tau) + \lambda \sum_{(\tau_i, \tau_j) \in \tau}K_p(\tau_i, \tau_j).$ Based on INViT-50-TSP, while maintaining the same experimental settings as in our paper, we trained the model with $\lambda$ set to 1, 0.1, and 0.01. A larger $\lambda$ indicates a higher importance of purity information in the loss function. The test results are as follows.

We can see that the test results for all three models are worse than those obtained with PUPO training. The degree of deterioration increases with increasing $\lambda$. Explicitly incorporating purity information into the loss function can lead to a mix of purity information and length information, which can lead to situations where the degree of purity improves but length deteriorates, making it impossible to distinguish well between them. This prevents the model from effectively perceiving purity information during training.

Reference:

[1] Reinforcement Learning: Theory and Algorithms. Alekh Agarwal, Nan Jiang, Sham M. Kakade, Wen Sun

Thank you very much for your time and effort in reviewing our work. We feel glad to know you acknowledge our work is a valuable contribution, find notion of purity is a sound underlying intuition, simple yet elegant, and believe this work would be insightful for the community.

We will address each of your concerns in detail as follow.

--Q1: Integration with REINFORCE.

Thank you for raising this concern. First, the proposed PUPO can be combined with any policy-based method. Due to the fundamental role of REINFORCE in NCO, we first combined PUPO with it and found that after PUPO-training, the model's cross-distribution and large-scale generalizabilities were significantly improved. Furthermore, we conducted experiments combining PUPO with imitation learning, and the results were further improved. (Please see our response to Reviewer Kzrh.)

In the field of neural routing problems, REINFORCE is a principal training method for neural constructive approaches, due to its ease of instance generation for the routing problems and its efficient sample utilization for reward computation. It is still adopted by the latest and most advanced works, such as [1-7].

Traditional reinforcement learning environments require CPU-based trajectory data generation, which is very time-consuming and slow. Therefore, advanced algorithms such as PPO are needed to improve sample utilization efficiency. However, in the neural routing problems, instance generation is very simple and easy to parallelize on GPUs, and reward computation can efficiently utilize samples. At the same time, [8] demonstrates that the complex hyper parameter tuning of PPO may slightly hurt performance. Also [8] shows that REINFORCE-trained model performs even better than PPO-trained on mDPP. Therefore, REINFORCE is a vital baseline for neural routing problem solvers training.

--Q2: Comparison against concorde.

Thank you for this suggestion. The test datasets used in this paper is from [6]. In the datasets with multi-distribution, TSP-100 is solved by Gurobi, an exact solver, so the solutions are guaranteed to be truly optimal. TSP-1000, TSP-5000 and TSP-10000 are solved by a SOTA heuristic algorithm LKH3, which can output near-optimal solutions.

To settle down your concerns, we use concorde to resolve three large datasets. We report the gap between LKH3 and concorde on TSP-1000 as follow, where negative gap means concorde performs better than LKH3.

It can be observed that LKH3 can reach the optima very nearly. However, due to unacceptable time consumption (one instance need more than 20 hours), concorde fails on the TSP-5000 and TSP-10000.

In the future, we will make the concorde-solved optimal values and solutions of each of 200 instances with four distributions in TSP-1000 all open source.

--Q3: Formal proof for purity phenomenon.

Thank you for raising this important concern. For the random Euclidean TSP, we theoretically prove that under certain conditions, the probability of an edge with purity order $k$ being at the optimal solution is bounded by $\frac{k!!}{4^{\frac{k}{2}-1}(k+1)!!}$. Notably, this upper bound decreases monotonically with increasing $k$. This means that edges of higher purity order have a lower probability of forming the optimal solution, consistent with our findings. Furthermore, the negative exponential form of this upper bound indirectly confirms the purity law we propose. If the paper is accepted, the theorem and detailed analysis and proof will be added to the paper.

THEOREM: For edge $e_{xy}$, assume that $K_p(e_{xy})=k$ is an positive even number, $z_i \in N_c(e_{xy})$ evenly distributed on both sides of $e_{xy}$, and the probability of each node being selected follows a uniform distribution. Then, the probability of $e_{xy}$ being at the optimal solution holds the following inequality:

$P[e_{xy} \in \tau_{OPT}] \leq \frac{k!!}{4^{\frac{k}{2}-1}(k+1)!!}.$

Let $\omega_{xy}$ denote the Hamilton Path formed by edge $e_{xy}$ and vertices $z_i \in N_c(e_{xy})$, $\tau_{\omega_{xy}}$ denote the Hamilton Circle including $\omega_{xy}$. Then, the following inequality holds:

$P[e_{xy} \in \tau_{OPT}] \leq P[$ edges in $\tau_{\omega_{xy}}$ do not cross $] \leq P[$edges in $\omega_{xy}$ do not cross $]$.

In the following, We will analyze the probability that no edges in $\omega_{xy}$ cross and obtain an upper bound. We prove the theorem by induction.

When $K_p(e_{xy})=k=2$, it can be proven by enumeration method.

Now, we first assume that when $K_p(e_{xy})=k$, the following equation holds:

P[ edges in $\omega_{xy}$ do not cross $] = \frac{k!!}{4^{k/2-1} (k+1)!!}$.

It means that among all $(k+1)!$ combinations of $\omega_{xy}$, there are $\frac{(k!!)^2}{4^{{k/2}-1}}$ situations where the edges of $\omega_{xy}$ do not cross.

When $K_p(e_{xy})=k+2$, first, there are totally $(k+3)!$ combinations of $\omega_{xy}$. Next, we consider the non-intersection structure based on case $k$. For every $\frac{(k!!)^2}{4^{{k/2}-1}}$ situations where the edges of $\omega_{xy}$ do not cross, consider one side of $e_{xy}$ first. If adding a point still leaves no intersection, there are $k/2+1$ possible addition positions. Considering both sides symmetrically, there are $(k/2+1)^2=\frac{(k+2)^2}{4}$ non-intersection cases. Therefore,

P[ edges in $\omega_{xy}$ do not cross $] = \frac{\frac{(k!!)^2}{4^{{k/2}-1}} \cdot \frac{(k+2)^2}{4}}{(k+3)!} = \frac{(k+2)!!}{4^{(k+2)/2-1} ((k+2)+1)!!}$.

Hence, we have finished the proof of $P[e_{xy} \in \tau_{OPT}] \leq \frac{k!!}{4^{\frac{k}{2}-1}(k+1)!!}$.

Reference:

[1] CaDA: Cross-Problem Routing Solver with Constraint-Aware Dual-Attention, ICML 2025

[2] SHIELD: Multi-task Multi-distribution Vehicle Routing Solver with Sparsity and Hierarchy, ICML 2025

[3] Rethinking Light Decoder-based Solvers for Vehicle Routing Problems, ICLR 2025

[4] PolyNet: Learning Diverse Solution Strategies for Neural Combinatorial Optimization, ICLR 2025

[5] UDC: A Unified Neural Divide-and-Conquer Framework for Large-Scale Combinatorial Optimization Problems, NeurIPS 2024

[6] INViT: A Generalizable Routing Problem Solver with Invariant Nested View Transformer, ICML 2024

[7] Towards Generalizable Neural Solvers for Vehicle Routing Problems via Ensemble with Transferrable Local Policy, IJCAI 2024

[8] RL4CO: an Extensive Reinforcement Learning for Combinatorial Optimization Benchmark

Dear Reviewer 8oG4,

We sincerely appreciate your thoughtful comments and suggestions again.

We would like to kindly follow up to check whether our response has addressed your concerns. If there are any remaining questions or suggestions, we would greatly appreciate the opportunity to clarify further.

If our responses have satisfactorily addressed your questions, we would be grateful if you would consider revisiting your rating.

Thank you once again for your time and constructive feedback.

Best regards,

The Authors of Submission 19649