Geometric Algorithms for Neural Combinatorial Optimization with Constraints

We would like to thank the reviewer for their detailed response.

While the authors make bold claims about "incorporating geometric algorithms in deep learning pipelines", all they do in the end is a standard Caratheodory decomposition.

The reviewer adds

but I would like to emphasize that this is not groundbreaking on a mathematical level and has been used in CombOpt algorithms for a long time.

• We believe that the statement of 'incorporating geometric algorithms in deep learning pipelines' is appropriate since the proposed method leverages constructive results from the literature on geometric algorithms to implement a differentiable pipeline for combinatorial optimization with neural nets. The theorem that forms the basis of our approach is in the GLS book 'Geometric algorithms and combinatorial optimization'. Our wording choice in the line quoted by the reviewer is a direct reference to the title of the book. We cite the book and mention the theorem explicitly in the paper (see line 192). We do not attempt to present this as some kind of outstanding mathematical achievement. The contribution lies in leveraging such algorithms for the purpose of neural comb opt. In that context, the constructive Caratheodory decomposition is able to effectively deal with both loss function design and rounding. We also provide extensive benchmarking to show the empirical benefits of the method (see also response to reviewer 7pF1). Our ablation also shows how this can even work using direct optimization with Adam.

• The reviewer acknowledges that such decompositions have been used in combinatorial algorithms for a long time. That is one of the strengths of this approach, since one can leverage results from the literature to adapt our method to other problems (e.g., scheduling, networking, etc.). To the best of our knowledge, the specific decompositions we describe have not been published but by no means constitute mathematical breakthroughs. That's because their existence essentially follows from the GLS result (theorem 4.1). They are valuable elements of our contribution, but they are not the central contribution of the paper. It is worth noting that the design of such efficient decompositions has been acknowledged as an interesting question in the literature (see for example section 4.3 in [1] or the introductory remarks in [2]). Introducing this problem as an element of machine learning algorithm design and making progress on it could be beneficial for both machine learning communities and communities working on algorithms. We certainly do not intend to mislead with our claims so we welcome any corrections of specific sentences that the reviewer believes would alleviate the issue.

Algorithm 1 is insufficiently described...

also

Theorem 4.1 contains many undefined terms and is not helpful in this formulation. I understand that it is an informal version...

• Here are are some clarifying comments that we will make sure to include in the final version of the paper:

  • By support of $\mathbf{x}$, we mean the nonzero coordinates of the vector.

  • The function $g$ is the (a.e.) differentiable function that computes the coefficient of the current iterate given the current iterate and a vertex of the polytope. We will make sure to introduce it earlier in the paper since this was also confusing to another reviewer. For example, the cardinality constraint section of the paper provides a concrete example for $g$ (see Equation 6).

  • Finding a vertex for the current iterate is problem dependent. In the most general case, according to the GLS theorem, we need to find a vertex of the smallest face of the polytope that contains the current iterate. This is further elaborated on in appendix B. For the problems we considered, a greedy algorithm can be sufficient. For instance, in the cardinality constraint we fetch the vertex corresponding to the the top-k coordinates. In some cases we may have to solve a linear program.

Why should it output a unique Caratheodory decomposition...

• The reviewer raises a valid point about the uniqueness of the decomposition. The uniqueness is required in order to ensure a.e. differentiability of the extension and we achieve uniqueness by designing a deterministic decomposition algorithm that maps a point from the polytope to a distribution. The details of such deterministic algorithm are problem specific and we provide several examples such as cardinality constraints and other matroids. For example, if we have a problem defined on graphs, we need to fix an ordering of the vertices of the input graph so as to consistently obtain a vertex and a coefficent corresponding to the current iterate. The decomposition then yields a deterministic mapping from the point in the polytope to the coefficients. In the case of the cardinality constraint, the set corresponding to the top-k coordinates of each iterate is chosen. Each iterate is a deterministic function of the previous iterate. Hence, the collection of sets (corresponding to top-k coordinates) and their coefficients are consistently obtained from the inputs. The ordering may only become relevant in case of tie-breaks (a rare occurrence in practice).

Shouldn't, by Caratheodory's theorem, the algorithm terminate after n+1 steps, with epsilon = 0, so why do you need this epsilon-test?

• The number of steps is always at most $n+1$ in theorem 4.1, which indeed has epsilon=0. Iterations should start from 1 and end at $n+1$ in the GLS theorem (Theorem 4.1). We apologize for the typo. Terms like strong optimization oracle are explained in the Appendix B but we will make sure to clear this up in the updated version.

• In appendix $C$ the role of the $\epsilon$ test is discussed in more detail. The main idea is to obtain a decomposition that contains potentially more than $n+1$ sets by rescaling the coefficient choice at every iteration. This allowed for better exploration (and convergence) during training and led to improved results. Furthermore, the $\epsilon$ test is motivated by the existence of approximate versions of the Caratheodory decomposition that do not depend on the ambient dimension of the space. We mention those in the last paragraph of the related work section of the main paper (references [3,12,47]). We consider this as a promising future direction for this line of work.

It is questionable whether the proposed pipeline is broadly applicable to difficult problems. The example problems in the paper (optimizing over a uniform matroid, partition matroid, spanning trees) are all easy problems from a complexity perspective.

• First, we want to acknowledge that depending on the problem, it may indeed be challenging to implement our approach. We mention this in the conclusion and in the extended discussion in Appendix B. We want to emphasize that cardinality constrained problems like maximum coverage or max-k-cut are known to be NP-Hard. However, the computational complexity of the problem is not sufficient to decide whether our method is applicable. As another example, there exist efficient decompositions for permutation polytopes and permutation polytopes can be used to solve hard (in terms of complexity) combinatorial problems like travelling salesperson. More generally, the viability of our approach primarily depends on whether one can find vertices and coefficients for each iteration of the decomposition efficiently. Hence, implementing our approach will be challenging for cases where linear optimization over the polytope is computationally hard. This is what is troublesome for independent set. We included the independent set example to highlight that even in cases where it is not obvious how to construct an efficient decomposition due to the potential intractability of the polytope, there are compromises that can be made to yield a viable method (i.e., a relaxation like FSTAB). FSTAB's corners are either integral vertices corresponding to independent sets or half-integral vertices. The decomposition will terminate but it may include some infeasible points (half-integral vertices). We may simply set the value of the objective to 0 for all non-integral points so that for an FSTAB extension the optimization will prioritize integral corners. Admittedly, this is not as mathematically clean as in other constraints but we wanted to bring up how one could deal with those harder cases.

I find it questionable to defer related work discussions to the appendix. If the work is relevant, put it in the main text. If not, leave it out.

• Given the space limitations, it is not viable to include almost 2 pages of related work in the main paper. We chose to prioritize the papers from the literature that we find are the most pertinent in the context of our contribution. The approach we take to those problems is not widely known, so our view was that the additional discussion around related work that we provide in the appendix could help interested readers contextualize our contribution and improve the overall accessibility of the paper.

As far as I understand, the authors have not submitted code for their experiments.

• We have submitted the code and the link can be found on page 17 of our corresponding NeurIPS checklist answer to this question. We will make sure to mention this in the updated main body of the paper too.

References:

1. Cunningham, William H. "On submodular function minimization." Combinatorica 5.3 (1985): 185-192.

2. Hoeksma, Ruben, Bodo Manthey, and Marc Uetz. "Efficient implementation of Carathéodory’s theorem for the single machine scheduling polytope." Discrete applied mathematics 215 (2016): 136-145.

Thank you for engaging with our rebuttal and responses!

the one thing that I am firm about and which prevents me from recommending acceptance is the insufficient presentation. I'd like to see clean presentation of Algorithm 1 and Thm. 4.1, and other places I commented on in my original review.

Could you please specify what is the exact issue with 4.1 or algorithm 1? Is there some specific problem with the answers we provided? To be clear, thm 4.1 in the paper is just a slightly abstracted version of the GLS theorem. This was done purely to save space. The full theorem from the book can be found in our appendix, together with the preliminaries required, but we reproduce the essential components here to facilitate discussion:

Let $\mathcal{P}\subseteq\mathbb{R}^n$ be a polyhedron and let $\varphi$ and $\nu$ be positive integers.

We say that $\mathcal{P}$ has facet‐complexity at most $\varphi$ if there exists a system of inequalities with rational coefficients that has solution set $\mathcal{P}$ and such that the encoding length of each inequality of the system is at most $\varphi$.

A well‐described polyhedron is a triple $(\mathcal{P}; n, \varphi)$ where $\mathcal{P}\subseteq\mathbb{R}^n$ is a polyhedron with facet‐complexity at most $\varphi$.

The strong optimization problem for a given rational vector $\mathbf{x}$ in $n$ dimensions and any well described $(\mathcal{P}; n, \varphi)$ is given by $\max{\mathbf{c}^\top \mathbf{x}}, \quad \mathbf{c} \in \mathcal{P}$.

A strong optimization oracle in our case is just an oracle for the above problem. The theorem is:

There exists a polynomial-time algorithm that for any well-described polyhedron $(\mathcal{P}; n, \varphi)$ given by a strong optimization oracle, for any rational vector $\mathbf{x}$, finds affinely independent vertices $\mathbf{1}_{S_1}, \mathbf{1} _{S_2} , \dots, \mathbf{1} _{S_k}$ and positive rational numbers $\lambda _1, \lambda _2, \dots, \lambda _k$ such that $\mathbf{x}= \sum _{t=1} ^k \lambda _t \mathbf{1} _ {S_t}$. Here $k\leq n+1$.

If the paper is accepted, we are given an additional page of content so including a few additional lines with definitions is easy. We want to emphasize that all of the above are contained in the standard textbook by GLS. It is difficult to provide a complete primer into polyhedral geometry within the page limit that covers all of the required background, since one could further inquire about the details of the encoding length, the precise notion of an oracle, etc.

Similarly for algorithm 1, could you please explain what the exact issue is?

If it is applicable to the matching problem, why haven't you mentioned this in the paper? Could you really solve, e.g., min-weight perfect matching in a general graph efficiently with your approach? Similarly for the approach you outline with respect to the TSP polytope: if this is a viable approach to improve existing TSP algorithms, why isn't it included in the paper?

There is no obstruction when it comes to applying our method to the matching problem or TSP. However, that is also the case for several other problems. In each case, we would have to provide some writeup (and ideally experiments) if we're going to present it. The paper runs over 40 pages already so we had to draw the line somewhere. For instance, suppose we decide to talk about TSP.

In that case, we could use the Birkhoff polytope, together with a Sinkhorn algorithm to obtain an interior point. Then we could use a fast matching oracle and results from [1] (for example) to perform the decomposition. The extension would then optimize over a distribution of tour lengths that is obtained via the distribution of permutations encoded by an interior point of the polytope. Obviously, while we can summarize this here in a few words, it would require a careful writeup in the appendix.

More generally, our approach makes more sense when no exact fast algorithm is available. If a problem is optimally solved by a fast greedy algorithm, then it's harder to justify using our approach.

I would be very curious to learn more about the polyhedral structure necessary to make this approach applicable.

The main sufficient condition is the existence of a polytime/fast oracle as seen in Theorem 4.1. With that it is possible to set up the algorithm. But as we discussed in the independent set example it is not necessary since one can work around it using a relaxation. How successful the approach is will depend on many implementation details such as the choice of model used (see also the k-cut ablation where the benefit of the model is shown).

We'll respond to the minor comments in a separate comment.

1. Dufossé, Fanny, et al. "Further notes on Birkhoff–von Neumann decomposition of doubly stochastic matrices." Linear Algebra and its Applications

We sincerely appreciate the reviewer’s thorough evaluation and constructive feedback. Below, we address the bulleted questions in the order they were presented.

• Thank you for pointing this out. You are absolutely right—the initialization should indeed be $\mathbf{x}_1\gets\mathbf{x}$, and we have corrected this accordingly. We also appreciate your suggestion regarding clarity: to improve readability, we will move the relevant definitions earlier in Section 4.1 so that Steps 3–5 are easier to follow when first encountered.

• The details of TTO are given in Lines 346 to 357. The caption of the Table has been updated to include the meaning of TTO (Test-Time Optimization). Table 1 is now renamed as Figure 1. The figures have been enlarged for readability; an enlarged version of the plots was already included in the appendix. In the last graph (12th), the overlap between Greedy and Gurobi is now explicitly noted.

• The MIP-based solution refers to the exact solution obtained by formulating the problem as a Mixed Integer Program (MIP) and solving it using Gurobi.

• A detailed discussion of how $\epsilon$ serves as a tolerance parameter to control the size of the decomposition is provided in Appendix Section C. The main idea is to 'undershoot' each step in the Caratheodory decomposition by rescaling the coefficient choice at every iteration in order to control the number of iterations and hence the number of sets that are used in the decomposition. This helped with convergence during training and led to improved results. Please see the discussion in Lines 320 to 329 around scaling factor and lines 1103 to 1110 in the Appendix for the choices of scaling factors. The scaling factors influence the number of terms generated in the decomposition—specifically larger scaling factor results in fewer terms in the decomposition. In the experiments $\epsilon$ is set to a very small value (e.g., $10^{-14}$) to account for numerical instability.

• In both Theorem 4.1 (due to Grötschel–Lovász–Schrijver) and our Theorem 4.2, we have $\epsilon=0$. As stated in Theorem 4.2, instantiating Algorithm 1 with Equations (5) and (6) yields an exact equality; that is, $\mathbf{x}=\sum_{t}p_{\mathbf{x}_t}(S_t)\mathbf{1} _{S_t}$.

We are thankful for the especially thorough review and their constructive comments. Below we respond to the bulleted questions in order.

• We tested two versions of the “short+Gurobi” approach on the random datasets, one with a time limit near our "medium" method, and for some settings, an "extended time" version, given the time close to our "long" method. When the instance size and the cardinality constraint are small (k = 10 on Random 500), "short+Gurobi" performs slightly better than other methods. For larger instances like Random 500 with k = 50 and Random 1000, our "medium" and "long" approaches show their advantages and achieved similar or better results, indicating that our local improvement technique effectively constrains the solution space and works better even compared to the well-optimized Gurobi Solver. In addition, the short+Gurobi method shows much better cost-efficiency compared to the pure Gurobi method, showing the potential of our method to serve as an efficient initial solution for structured solvers, significantly reducing overall computation time.

Dataset: Random500

* Gurobi finishes search in less than 5 seconds, so extended time version is not applicable

Dataset: Random1000

• In general, Carathéodory’s theorem applies to any convex hull. However, the specifics of obtaining an almost-everywhere (a.e.) differentiable extension can vary when moving beyond the Boolean domain and may require domain-specific techniques rather than direct generalization. In particular, Steps 3, 5 and 6 of our Algorithm 1 are tailored to the $\{0,1\}$ setting and may require modification in other domains. These steps are problem-specific, and obtaining analogous steps in non-Boolean settings would be a worthwhile contribution on its own. For example, in combinatorial problems over permutations—where we are no longer dealing with $\{0,1\}$ assignments—achieving an a.e. differentiable extension typically relies on the Birkhoff decomposition. This approach requires additional problem-specific algorithmic tools. Recent work by Nerem et al. [1] explores this direction in more depth.

• Thank you for the constructive feedback on the presentation—we sincerely appreciate it. We will incorporate your suggestions in the revised version of the paper.

  • Clarification on the meaning of each "X" in the plots: In each plot, the x-axis represents runtime (where smaller values are better), and the y-axis represents the objective value (where larger values are better). Each colored cross ("X") corresponds to a method identified in the legend and indicates that method’s performance in terms of both runtime and objective. The relative positions of the crosses illustrate the trade-off between these two metrics. Across the various settings, our methods consistently appear on or near the Pareto front, demonstrating favorable performance in both dimensions.
  • We will change both the marker shapes and colors to improve readability.
  • We agree that standardizing the axis ranges can improve visual clarity. However, since we experiment with different values of $k$ (in the cardinality constraint), the optimal objective value naturally increases with larger $k$, which leads to differing y-axis scales across plots. This variation reflects meaningful differences in the problem instances rather than plotting inconsistencies. That said, we appreciate the point about the Random (240s) baseline. Given its poor performance and the distortion it introduces in the axis scaling, we will move it to the appendix and add a note in the main text to guide the reader. We also note that larger versions of all plots are already available in the appendix for more detailed comparison.
  • We will clean up Tables 9–11 and consider splitting the tables or reducing precision to improve readability.
  • The reason the reference numbers in the Appendix differ from those in the main paper is that we include additional related work and citations in the Appendix (Section A) that are not referenced in the main text. As a result, the numbering continues beyond what appears in the main paper. There are a few inconsistencies that we have fixed. Thanks for pointing this out.

References:

1. Robert R Nerem, Zhishang Luo, Akbar Rafiey, and Yusu Wang. Differentiable extensions with rounding guarantees for combinatorial optimization over permutations. arXiv preprint arXiv:2411.10707, 2024

We would like to thank the reviewer for their thoughtful and constructive feedback on our work.

• The reviewer makes a point about the applicability of the method to problems where feasible solutions are difficult to obtain. Indeed, not being able to efficiently construct feasible solutions poses a challenge and our proposed method works best when obtaining feasible solutions is computationally tractable. This allows us to construct a decomposition by fetching feasible sets at each iteration. On the other hand, there are plausible ways in which this approach can be naturally extended to the setting where feasible solutions are hard to construct. For instance, one could define a polytope of 'efficiently computable infeasible solutions' and optimize over them by scoring the amount of constraint violation. This could be used as warm-starts that are passed to a discrete search algorithm that then attempts to construct a feasible solution. This is an interesting future direction but it introduces an additional layer of complexity in terms of benchmarking and implementation so we decided to prioritize simpler problems.

The reviewer also asked for a summary of the results in the appendix, particularly for problems other than max coverage.

• To summarize the results in the appendix, our method achieves strong empirical performance on the max coverage problem and is competitive with the greedy approach on the max k-cut problem. The appendix also includes several ablation studies that demonstrate the effectiveness of our medium and long methods, as well as the quality of candidate sets returned by the trained model. Additionally, we highlight our focus on the cardinality constraint and discuss the shortcomings of the best-known neural baseline.

• Below, we expand on these points.

  • While the appendix presents extensive experiments on the max coverage problem—using various cardinality constraints and datasets—as well as results on the max k-cut problem, we are also in the process of obtaining additional results for other problems. Below we provide preliminary results for another fundamental combinatorial problem: the Facility Location problem. Here we present the performance of our short method (with no Test-Time Optimization) against the Greedy, an RL baseline, and the CardNN approach without Test-Time Optimization (TTO). We consider two cardinality constraints cases that are k=5 and k=20 and report the objective function values. Note that Facility Location is a minimization problem and smaller objective values are better. The experiments are run on synthetic data where each instance contains 500 points.

  Method	k=5	k=20
  Ours short	16.1775	4.1373
  CardNN-S-NoTTO-NoKmeans	86.2816	17.6436
  Greedy	17.3714	2.9644
  RL	35.9455	6.8085

  • In Appendix G.3, we run an ablation study to demonstrate that in the medium and long versions of our approach, the strategic selection of nodes in our candidate pool provides meaningful benefits. The results in Table 2 demonstrate that our strategic candidate selection provides consistent improvements over random selection, validating that our local improvement procedure benefits specifically from the candidate pool construction via our decomposition, rather than simply from having access to additional nodes for local improvement.

  • Note that we heavily emphasized cardinality constraints in the experiments because they contain practically important problems like max-coverage. Despite the simplicity of the constraint, we found that existing solutions in the neural combinatorial optimization literature did not provide sufficiently simple and/or effective ways of tackling the constraint while achieving strong results. The ablation study we performed in page 31 (Appendix G.4) highlights this issue with the best known neural baseline. Through our projection-decomposition-extension scheme, we provide an intuitive and empirically compelling approach to solving the problem in a differentiable way, and our experimental results support that.

It would be better to have a numerical representation of those figures. In addition, the result tables in appendix are how to read as well.

• The reviewer also raises a valid point about the presentation of the results which we completely acknowledge. We conducted a comprehensive set of experiments on our method but also on the competing baselines to ensure a fair evaluation. There are several components to those comparisons. Different datasets, different parameters for the constraint, as well as the difference between baselines that involve learning and those that don't. While preparing the paper, it was not clear to us what is the best way to present these results, because the full table (as can be seen in the appendix G.5) is also quite dense and hard to read. We are going to do our best to reformat the results to highlight some of the key takeaways according to the reviewer's comment, but we welcome any further recommendations on how we could present the results.

• We will change both the marker shapes and colors to improve readability. We will clean up Tables 9–11 and consider splitting the tables or reducing precision to improve readability. We also note Appendix G.5 contains larger versions of all plots for more detailed comparison as well as raw numbers for experiments in tables 9 to 11.

Could you comment on how well this approach scale to larger instances? How computationally heavy is differentiating through the decomposition steps?

• Regarding computational efficiency, the detailed cost depends on the problem and the cost of the decomposition, which in turn depends on the algorithm used to obtain feasible sets and the ability to compute intersections with the boundary of the polytope. Given a point in the polytope, the decomposition applies iteratively a transformation to the point until the termination criterion is reached. For the examples we discussed, that is roughly $O(n)$ (in the worst case) iterations given by Equation (3). For the backprop step, the gradients are passing through the coefficients of the decomposition, where we have one coefficient per iterate. Therefore, the cost of gradient computation for each coefficient depends linearly on the number of iterations it took to compute the coefficient because the chain rule for a given coefficient will have to go through all the iterates before it. Crucially, the most costly part is actually the decomposition itself and the backward pass does not increase the overall time significantly. As a general rule, the closer the neural net prediction is to a feasible solution, the fewer the decomposition steps and hence the method becomes more efficient. How close the predictions are to feasible depends on practical implementation, how the architecture aligns with the problem, etc. We were able to run our experiments containing instances with several thousands of nodes on a single 48GB GPU. To make the point about the cost of backprop clearer, we provide some time measurements in the coverage problem.

• For a given instance size, the decomposition will run at most $n-k$ times which explains why for smaller values of $k$ the cost is larger. Overall, we see that the cost of backprop is not significantly increased as we scale up the instances.

Why did you choose to present ablation study results on max-cut which is different from the coverage problem?

• Finally, the reason for choosing max k-cut for the ablation was to highlight another problem where our approach can work. It is also worth noting that the same ablation shows that our decomposition-based extension also works with direct optimization on the input space (i.e., without a neural network). Even though the constructed function is non-convex, one can just directly run Adam on a random starting point in the hypercube to find a solution for the problem. If the paper is accepted, we are planning to include even more results along those lines to further bolster our claims.