Self-Supervised Combinatorial Optimization with Constraints via Frank-Wolfe

# Self-Supervised Combinatorial Optimization with Constraints via Frank–Wolfe
Source: [https://arxiv.org/html/2609.25728](https://arxiv.org/html/2609.25728)
Yifei Xu11footnotemark:1Affiliation:NYUEmail:[yx3590@nyu\.edu](mailto:)Nikolaos KaraliasAffiliation:MITEmail:[stalence@mit\.edu](mailto:)

###### Abstract

Self\-supervised learning for combinatorial optimization has emerged as a promising paradigm for solving discrete optimization problems with neural networks, but a central challenge remains: handling hard combinatorial constraints within continuous, gradient\-based training\. Continuously extending combinatorial objectives to convex domains is a powerful technique, yet existing approaches often require projection steps that constrain neural network outputs to lie inside the feasible polytope and rely on ad\-hoc and problem\-specific constructions\. We propose a general framework in which the neural network is allowed to predict arbitrary continuous vectors that could potentially lie outside of the feasible polytope\. These predictions are then approximated by sparse convex combinations of feasible solutions using a geometric decomposition algorithm based on Frank–Wolfe methods and approximate Carathéodory results\. This decomposition induces an a\.e\.\-differentiable, self\-supervised loss defined as the expected value of the discrete objective\. The same procedure provides an automatic rounding guarantee at inference time\. We demonstrate strong empirical performance across multiple combinatorial problems, including the Quadratic Assignment Problem, Maximum Coverage, and the Traveling Salesperson Problem\.

## 1Introduction

Combinatorial optimization \(CO\) forms a central pillar of optimization theory and practice, encompassing a broad class of problems in which one seeks to optimize an objective over a discrete, often exponentially large, feasible set\. Such problems arise ubiquitously across science and engineering, from routing and resource allocation to learning and inference\. Canonical examples include the Traveling Salesperson Problem \(TSP\), Maximum Coverage \(MC\), and Quadratic Assignment Problem \(QAP\), among many others\. What unifies these problems is not only their computational hardness, but also the rich algorithmic structure induced by discrete constraints that describe large\-scale discrete combinatorial spaces of configurations\.

Despite this shared structure, successful approaches to CO have historically been case\-specific, often reflecting deep problem\-dependent insights; "stroke\-of\-genius" heuristics\. A common algorithmic template, especially in approximation algorithms, relaxes the discrete problem into a continuous or convex surrogate, solves for a fractional solution, and then applies a bespoke rounding or improvement procedure to recover feasibility\. While this relax\-and\-round paradigm has led to remarkable practical results and theoretical guarantees, adapting this approach to gradient\-based data\-driven settings can be challenging\. This limitation and challenge is especially pronounced when problem instances are drawn from structured distributions, increasingly common in real world applications\. Self\-supervised learning \(SSL\) approaches for combinatorial optimization offer an appealing alternative: unlike static approximation algorithms and heuristics, they can leverage latent patterns in data, train efficiently without requiring large amounts of labeled solutions, and incorporate algorithmic priors that facilitate learning and generalization\.

In this work, we bridge the gap between discrete optimization and deep learning by building on work that used continuous extensions of discrete functions as losses for neural CO\[[24](https://arxiv.org/html/2609.25728#bib.bib16),[23](https://arxiv.org/html/2609.25728#bib.bib11),[50](https://arxiv.org/html/2609.25728#bib.bib21)\]\. The central idea in this line of work is to use extensions to embed discrete constraints and objectives directly into a learning pipeline, allowing learning and optimization to proceed end\-to\-end\. By smoothing the discrete landscape, these approaches enable an unsupervised learning paradigm where models learn to exploit patterns and discover high\-quality solutions without requiring expensive labels\.

More specifically, we propose a generic learning\-based pipeline that integrates hard combinatorial constraints directly into training and inference\. The neural network is allowed to output an arbitrary continuous vector, without any feasibility requirement\. This output is then passed through a geometric decomposition algorithm inspired by approximate Carathéodory and Frank–Wolfe methods\[[4](https://arxiv.org/html/2609.25728#bib.bib12),[46](https://arxiv.org/html/2609.25728#bib.bib1),[12](https://arxiv.org/html/2609.25728#bib.bib5)\], which produces a sparse convex combination of feasible solutions\. The resulting decomposition may differ substantially from the network output but is always expressed entirely in terms of feasible points\. This effectively shifts the burden of constraint enforcement from tuning penalty terms to capturing constraints algorithmically through our decomposition algorithm\. The decomposition induces a \(a\.e\. differentiable\) distribution over feasible solutions, and the expected value of the discrete objective under this distribution is used directly as a self\-supervised training loss\. In addition, the discrepancy between the network output and its decomposition can be quantified and included as an auxiliary regularization term, with an associated weight\. While not required by the framework, we find empirically that including this term improves optimization stability and solution quality\. Gradients propagate through this objective using automatic differentiation\. At inference time, the same procedure yields a small set of feasible candidate solutions from which the best one is selected\. Figure[1](https://arxiv.org/html/2609.25728#S1.F1)summarizes the pipeline\.

Our pipeline is projection\-free\. The decomposition in\[[23](https://arxiv.org/html/2609.25728#bib.bib11)\]requires the neural prediction to lie inside the feasible polytope, with feasibility ensured separately through problem\-specific constructions\. For example, it uses distinct projections for the hypersimplex \(Proposition 4\.4\) and the spanning\-tree polytope \(Theorem E\.7\)\. Our method removes the need for projection and accepts arbitrary ambient\-space outputs and uses a Frank–Wolfe decomposition to construct a sparse convex combination of feasible vertices\.

Moreover, our decomposition has an explicit run time certificate and is more efficient in terms of the number of Linear Maximization Oracle \(LMO\) calls\. It requires exactly one exact or approximate LMO call per iteration, soTTiterations use exactlyTTLMO calls\. In contrast, each iteration of the generic GLS \(Groetschel/Lovasz/Schrijver\) decomposition in\[[23](https://arxiv.org/html/2609.25728#bib.bib11)\]invokes minimal\-face and ray\-boundary optimization routines that themselves require solving convex optimization problems using the oracle; see the proof of Theorem 4\.1\. Consequently, a single GLS iteration may require many LMO calls\. For example, a standard ellipsoid\-based implementation requires approximatelyO⁡\(n2​L\)O\(n^\{2\}L\)oracle calls per iteration andO⁡\(n3​L\)O\(n^\{3\}L\)calls overall, whereLLis a precision parameter\.

Although specialized GLS decompositions can be more efficient for particular polytopes, deriving them requires nontrivial, problem\-specific techniques\. Given an efficient exact or approximate LMO, our procedure applies unchanged across different polytopes, with only the LMO changing\. It thus provides a unified oracle\-based method for hypersimplex, Birkhoff, spanning\-tree, and matching polytopes, subsuming the specialized settings considered in\[[23](https://arxiv.org/html/2609.25728#bib.bib11),[50](https://arxiv.org/html/2609.25728#bib.bib21)\]\.

![Refer to caption](https://arxiv.org/html/2609.25728v1/Figures/overall-FW.png)Figure 1:Overview of our framework\. During training, the model with parametersθ\\thetaoutputs a point𝐱θ\\mathbf\{x\}\_\{\\theta\}, which may lie outside𝒫\\mathcal\{P\}, the convex hull of feasible solutions\. Our Frank–Wolfe decomposition layer constructs a feasible proxy𝐮\\mathbf\{u\}as a convex combination of at mostTTfeasible vertices𝐯t\\mathbf\{v\}\_\{t\}, with weightsαt​\(𝐱θ\)\\alpha\_\{t\}\(\\mathbf\{x\}\_\{\\theta\}\)\. The self\-supervised lossℒSSL\\mathcal\{L\}\_\{\\mathrm\{SSL\}\}is computed from the expected discrete objectiveffand a reconstruction term weighted byλ≥0\\lambda\\geq 0\. The weightsαt​\(𝐱θ\)\\alpha\_\{t\}\(\\mathbf\{x\}\_\{\\theta\}\)are a\.e\. differentiable w\.r\.t𝐱θ\\mathbf\{x\}\_\{\\theta\}, allowing gradients from the loss to propagate through the decomposition layer back to the neural network\. During inference \(green box\), the same procedure generates feasible candidates, and the best solution𝐯∗\\mathbf\{v\}^\{\*\}is selected\. Solid arrows indicate the forward pass, and the dashed arrow indicates backpropagation\.#### Contributions\.

We make the following contributions:

- •A projection\-free, self\-supervised framework for CO\.We introduce a unified learning framework that enables neural networks to optimize combinatorial objectives under hard discrete constraints by operating over convex polytopes of feasible solutions\. Unlike prior geometric\-extension methods, our framework does*not*require network outputs to lie in the polytope during training and provides a common formulation across a broad class of combinatorial structures\.
- •An oracle\-based extension layer using Frank–Wolfe decomposition\.To realize this framework, we propose a Frank–Wolfe–style extension layer that maps an arbitrary continuous prediction to a sparse convex combination of feasible solutions, yielding a distribution over feasible discrete solutions\. Given an efficient exact or approximate linear maximization oracle \(LMO\), the same procedure applies across different polytopes; only the implementation of the LMO changes\. It requires one LMO call per decomposition iteration and avoids problem\-specific projections and specialized decomposition algorithms\.
- •A differentiable, self\-supervised loss with built\-in rounding guarantees\.The induced convex decomposition allows us to define a differentiable loss as the expected discrete objective, with approximation error explicitly incorporated\. At inference time, selecting the best solution in the support yields a feasible solution whose objective value is no worse than the expectation\. This provides a direct rounding guarantee and unifies learning and solution recovery within a single pipeline\.
- •Empirical validation across diverse combinatorial domains\.We demonstrate the effectiveness of our framework on a range of classical CO problems—including MC, TSP, and QAP—showing strong empirical performance and broad applicability across diverse constraint families\.

## 2Related Work

Extensions and optimization\.A long line of work in CO views discrete problems through the lens of continuous optimization by embedding feasible solutions into a convex polytope and optimizing over this relaxation\. This perspective is deeply rooted in algorithm design and polyhedral geometry, where constructing convex or concave extensions with favorable optimization properties has been a central theme\[[43](https://arxiv.org/html/2609.25728#bib.bib26),[13](https://arxiv.org/html/2609.25728#bib.bib27),[48](https://arxiv.org/html/2609.25728#bib.bib24),[60](https://arxiv.org/html/2609.25728#bib.bib25)\]\. In particular, the convex closure provides the tightest convex extension of a discrete set function and underlies many relaxation\-based methods\[[16](https://arxiv.org/html/2609.25728#bib.bib28),[60](https://arxiv.org/html/2609.25728#bib.bib25)\]\. Classical solvers and approximation algorithms routinely exploit this framework: prominent examples include LP and SDP relaxations for problems such as Max Cut and TSP, with solvers like Concorde and Gurobi relying on large\-scale LP and cutting\-plane techniques\[[19](https://arxiv.org/html/2609.25728#bib.bib29),[67](https://arxiv.org/html/2609.25728#bib.bib30),[14](https://arxiv.org/html/2609.25728#bib.bib31),[2](https://arxiv.org/html/2609.25728#bib.bib32)\]\. While these approaches yield strong theoretical guarantees, they typically rely on problem\-specific, rounding procedures that are non\-differentiable and difficult to integrate into modern, data\-driven learning pipelines\.

Neural CO\.Recent work has explored unsupervised and self\-supervised learning approaches for CO that optimize continuous surrogates of discrete objectives, typically induced through probabilistic relaxations and expected\-value formulations\[[1](https://arxiv.org/html/2609.25728#bib.bib13),[22](https://arxiv.org/html/2609.25728#bib.bib14),[61](https://arxiv.org/html/2609.25728#bib.bib15),[63](https://arxiv.org/html/2609.25728#bib.bib17)\]\. Compared to reinforcement learning or supervised methods\[[28](https://arxiv.org/html/2609.25728#bib.bib20),[31](https://arxiv.org/html/2609.25728#bib.bib34),[62](https://arxiv.org/html/2609.25728#bib.bib33)\], these approaches avoid labeled data and often exhibit more stable training\. Unsupervised CO pipelines involve several crucial components, including neural network architecture\[[45](https://arxiv.org/html/2609.25728#bib.bib35),[62](https://arxiv.org/html/2609.25728#bib.bib33),[54](https://arxiv.org/html/2609.25728#bib.bib36)\]and the role of input features for well\-known classes of models\[[40](https://arxiv.org/html/2609.25728#bib.bib37),[27](https://arxiv.org/html/2609.25728#bib.bib38)\]\. In this work, our focus is on loss function design and rounding\. A common strategy in self\-supervised CO is to parameterize a distribution over discrete solutions and use the expected value of the discrete objective as a differentiable training loss\[[22](https://arxiv.org/html/2609.25728#bib.bib14),[63](https://arxiv.org/html/2609.25728#bib.bib17),[7](https://arxiv.org/html/2609.25728#bib.bib39)\], enabling smoother optimization but raising challenges related to constraint enforcement\. A major consideration in these approaches is enforcing constraints on the outputs of neural networks\. Techniques in the literature typically involve projecting the network output onto the feasible set\. This has been explored for cardinality constraints using ideas from optimal transport\[[66](https://arxiv.org/html/2609.25728#bib.bib41)\], and for general linear constraints using gradient\-based methods\[[15](https://arxiv.org/html/2609.25728#bib.bib42),[69](https://arxiv.org/html/2609.25728#bib.bib40)\]\. This also includes Sinkhorn\-style projection methods and their extensions\[[56](https://arxiv.org/html/2609.25728#bib.bib49),[64](https://arxiv.org/html/2609.25728#bib.bib50)\]\. In the continuous \(and potentially non\-convex\) optimization setting, several projection techniques have also been developed\[[37](https://arxiv.org/html/2609.25728#bib.bib45),[38](https://arxiv.org/html/2609.25728#bib.bib46),[39](https://arxiv.org/html/2609.25728#bib.bib47),[36](https://arxiv.org/html/2609.25728#bib.bib48)\]\. Other approaches add penalty terms or improve training dynamics through annealing or physics\-inspired formulations\[[57](https://arxiv.org/html/2609.25728#bib.bib43),[55](https://arxiv.org/html/2609.25728#bib.bib44),[47](https://arxiv.org/html/2609.25728#bib.bib18)\]\. Our method differs from these approaches by incorporating feasibility directly through geometric decomposition\. We define a distribution supported entirely on feasible solutions, allowing constraints to be built into the distribution itself\. Thus, the training loss focuses solely on the discrete objective while remaining a\.e\. differentiable\. A closely related line of work develops neural CO methods based on geometric relaxations of feasibility polytopes\.[Karalias et al\. \[24\]](https://arxiv.org/html/2609.25728#bib.bib16),[Karalias et al\. \[23\]](https://arxiv.org/html/2609.25728#bib.bib11)propose differentiable extensions for set\-function optimization under various constraints, while[Nerem et al\. \[50\]](https://arxiv.org/html/2609.25728#bib.bib21)study CO problems over permutations\. These approaches typically require the neural network output to lie within the feasibility polytope\. In contrast, our approach allows arbitrary continuous network outputs and applies a general\-purpose geometric decomposition across a broad class of combinatorial constraints\.

Approximate Carathéodory and Frank\-Wolfe\.Our approach builds on results from convex geometry and projection\-free optimization\. A central component is the Frank–Wolfe algorithm\[[17](https://arxiv.org/html/2609.25728#bib.bib6),[9](https://arxiv.org/html/2609.25728#bib.bib2)\], a first\-order method for constrained optimization\. Several variants with improved convergence have been studied\[[21](https://arxiv.org/html/2609.25728#bib.bib3),[35](https://arxiv.org/html/2609.25728#bib.bib7)\], and Frank–Wolfe has been widely applied to structured and combinatorial polytopes\[[18](https://arxiv.org/html/2609.25728#bib.bib9),[34](https://arxiv.org/html/2609.25728#bib.bib8),[32](https://arxiv.org/html/2609.25728#bib.bib10)\]\. Approximate variants of Frank–Wolfe with provable guarantees have also been proposed\[[42](https://arxiv.org/html/2609.25728#bib.bib51),[70](https://arxiv.org/html/2609.25728#bib.bib4)\]\. From a geometric perspective, the Carathéodory theorem\[[10](https://arxiv.org/html/2609.25728#bib.bib22)\]states that any point in a polytopeP⊆ℝdP\\subseteq\\mathbb\{R\}^\{d\}can be expressed as a convex combination of at mostd\+1d\+1corners of the polytope\. Approximate versions provide dimension\-independent bounds on sparse approximations\[[4](https://arxiv.org/html/2609.25728#bib.bib12),[46](https://arxiv.org/html/2609.25728#bib.bib1)\]\. Obtaining such sparse convex combinations of polytope vertices is essential for our framework, and Frank–Wolfe\-type algorithms provide a natural algorithmic mechanism to achieve this\[[46](https://arxiv.org/html/2609.25728#bib.bib1),[12](https://arxiv.org/html/2609.25728#bib.bib5)\]\.

## 3Problem Formulation and Learning Setup

We consider CO problems defined over a discrete feasible set\. Each problem instance is specified by a collection of variablesV=\{v1,…,vn\}V=\\\{v\_\{1\},\\dots,v\_\{n\}\\\}, each assigned a value from a discrete domainD=\{1,…,d\}D=\\\{1,\\dots,d\\\}, a finite set of feasible solutions𝒞⊆Dn\\mathcal\{C\}\\subseteq D^\{n\}, and a real\-valued objective functionf:𝒞→ℝf:\\mathcal\{C\}\\to\\mathbb\{R\}\. The goal is to find an optimal feasible solution𝐱∗∈𝒞\\mathbf\{x\}^\{\*\}\\in\\mathcal\{C\}that optimizes \(either minimizes or maximizes\) the objective function:

max⁡f⁡\(𝐱\),s\.t\.𝐱∈𝒞\.\\displaystyle\\max f\(\\mathbf\{x\}\),\\quad\\text\{s\.t\.\}\\quad\\mathbf\{x\}\\in\\mathcal\{C\}\.\(1\)This abstract formulation accommodates a broad class of CO problems, including those defined over binary vectors, general discrete assignments, and permutations\.

Many examples in this paper focus on problems with Boolean decision variables, where solutions can be represented by indicator vectors encoding subset or structure selection under constraints\. Two representative running examples are Maximum Coverage, which selects a subset of elements subject to a budget constraint, and the TSP, where binary variables indicate whether edges are included in a Hamiltonian tour\. Our framework also applies to CO problems over permutations, e\.g\., the QAP, and to non\-Boolean discrete domains, e\.g\.,kk\-submodular maximization\.

A standard geometric perspective on CO is to consider the convex hull of all feasible solutions,

𝒫=conv⁡\(𝒞\)\.\\mathcal\{P\}=\\mathrm\{conv\}\(\\mathcal\{C\}\)\.For CO problems with Boolean domain, we have𝒫=conv⁡\(𝒞\)⊆\[0,1\]n\\mathcal\{P\}=\\mathrm\{conv\}\(\\mathcal\{C\}\)\\subseteq\[0,1\]^\{n\}\. As an example, when𝒞\\mathcal\{C\}consists of all subsets of sizekk,𝒫\\mathcal\{P\}is the associated cardinality polytope, a\.k\.a hypersimplex\. EachS∈𝒞S\\in\\mathcal\{C\}is represented by its indicator vector𝟏S∈\{0,1\}n\\mathbf\{1\}\_\{S\}\\in\\\{0,1\\\}^\{n\}, where\(𝟏S\)i=1\(\\mathbf\{1\}\_\{S\}\)\_\{i\}=1if elementiibelongs toSSand\(𝟏S\)i=0\(\\mathbf\{1\}\_\{S\}\)\_\{i\}=0otherwise\. The polytope associated with𝒞\\mathcal\{C\}is then𝒫:=conv⁡\{𝟏S:S∈𝒞\}⊆\[0,1\]n,\\mathcal\{P\}:=\\operatorname\{conv\}\\\{\\mathbf\{1\}\_\{S\}:S\\in\\mathcal\{C\}\\\}\\subseteq\[0,1\]^\{n\},whose vertices, i\.e\., extreme pointsext⁡\(𝒫\)\\operatorname\{ext\}\(\\mathcal\{P\}\), coincide with the indicator vectors of subsets of sizekk\.

Note that although𝒫\\mathcal\{P\}may have an exponentially large facet description, many important combinatorial polytopes admit efficient linear optimization via polynomial\-time separation or optimization oracles\.

### 3\.1Self\-Supervised Learning for CO

In the self\-supervised learning setting, each problem instanceℐ\\mathcal\{I\}is described by input features𝐙ℐ\\mathbf\{Z\}\_\{\\mathcal\{I\}\}possibly together with instance\-specific structure\. These features are processed by a neural networkNNθ\\mathrm\{NN\}\_\{\\theta\}, which outputs a representation𝐱θ\\mathbf\{x\}\_\{\\theta\}\. Importantly, this output is not required to lie in the convex polytope of feasible solutions\.

Rather than training the network to imitate optimal solutions, learning is guided by a loss function defined on an a\.e\. differentiable extension of the combinatorial objective\. This extension is constructed using the geometry of the polytope𝒫\\mathcal\{P\}, allowing the objective to be meaningfully evaluated on fractional solutions and optimized using gradient\-based methods\.

The network is trained over a collection of instances without access to ground\-truth optimal solutions, making the design of the loss function central\. At inference time, the continuous output𝐱\\mathbf\{x\}is converted into a feasible discrete solution using a principled, geometry\-based procedure, rather than heuristic rounding\. Moreover, because no labeled data is required, the same objective can be optimized directly on test instances, allowing additional computation at test time to further improve solution quality\.

## 4Proposed Method

We now describe the overall pipeline\. Consider a polytope𝒫=conv⁡\(𝒞\)\\mathcal\{P\}=\\mathrm\{conv\}\(\\mathcal\{C\}\), where𝒞⊂ℝn\\mathcal\{C\}\\subset\\mathbb\{R\}^\{n\}denotes the set of feasible solutions for a CO problem\. Letext⁡\(𝒫\)⊆𝒞\\operatorname\{ext\}\(\\mathcal\{P\}\)\\subseteq\\mathcal\{C\}denote the set of extreme points i\.e\., corners, of𝒫\\mathcal\{P\}which correspond to feasible combinatorial solutions\.

Given a target point𝐱∈ℝn\\mathbf\{x\}\\in\\mathbb\{R\}^\{n\}, which can be thought of as the output of a neural network, our goal is to find a point𝐮∈𝒫\\mathbf\{u\}\\in\\mathcal\{P\}that serves as a good proxy for𝐱\\mathbf\{x\}and satisfies desirable properties\. In particular, we want𝐮\\mathbf\{u\}to be a*sparse*convex combination ofTTpoints𝐯0,…,𝐯T−1∈ext⁡\(𝒫\)\\mathbf\{v\}\_\{0\},\\dots,\\mathbf\{v\}\_\{T\-1\}\\in\\operatorname\{ext\}\(\\mathcal\{P\}\)\. That is, for a smallTT, we have

𝐮⁡\(𝐱\)=∑t=0T−1αt​\(𝐱\)​𝐯t∈𝒫with\{𝐯t∈ext⁡\(𝒫\),∑t=0T−1αt\(𝐱\)=1;αt≥0\.\\displaystyle\\mathbf\{u\}\(\\mathbf\{x\}\)=\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\(\\mathbf\{x\}\)\\,\\mathbf\{v\}\_\{t\}\\in\\mathcal\{P\}\\quad\\text\{with\}\\quad\\begin\{cases\}\\mathbf\{v\}\_\{t\}\\in\\operatorname\{ext\}\(\\mathcal\{P\}\),\\\\ \\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\(\\mathbf\{x\}\)=1;~\\alpha\_\{t\}\\geq 0\.\\end\{cases\}\(2\)
The coefficientsαt​\(𝐱\)\\alpha\_\{t\}\(\\mathbf\{x\}\)naturally define a distribution𝒟⁡\(𝐱\)\\mathcal\{D\}\(\\mathbf\{x\}\)over feasible solutions\. We define the FW\-induced objective

F⁡\(𝐱\)=∑t=0T−1αt​\(𝐱\)​f​\(𝐯t\)=𝔼𝐯∼𝒟⁡\(𝐱\)​\[f⁡\(𝐯\)\]\.\\displaystyle F\(\\mathbf\{x\}\)=\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\(\\mathbf\{x\}\)\\,f\(\\mathbf\{v\}\_\{t\}\)=\\mathbb\{E\}\_\{\\mathbf\{v\}\\sim\\mathcal\{D\}\(\\mathbf\{x\}\)\}\[f\(\\mathbf\{v\}\)\]\.\(3\)
Training\.In the learning pipeline, let𝐱θ\\mathbf\{x\}\_\{\\theta\}denote the prediction of a neural network with parametersθ\\theta\. For a maximization problem \([1](https://arxiv.org/html/2609.25728#S3.E1)\), we train the network by minimizing the self\-supervised loss

ℒSSL​\(θ\):=−F⁡\(𝐱θ\)\+λ​‖𝐱θ−𝐮⁡\(𝐱θ\)‖22,λ≥0\.\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{SSL\}\}\(\\theta\):=\-F\(\\mathbf\{x\}\_\{\\theta\}\)\+\\lambda\\\|\\mathbf\{x\}\_\{\\theta\}\-\\mathbf\{u\}\(\\mathbf\{x\}\_\{\\theta\}\)\\\|\_\{2\}^\{2\}\\,,\\qquad\\lambda\\geq 0\.\(4\)
In order to optimize the loss function using gradient\-based methods, the weightsαt​\(𝐱θ\)\\alpha\_\{t\}\(\\mathbf\{x\}\_\{\\theta\}\)must be differentiable functions of𝐱θ\\mathbf\{x\}\_\{\\theta\}\. Under this condition, derivatives of the loss with respect to the network parameters can be computed via the chain rule:∂F∂θ=∂F∂𝜶⁡\(𝐱θ\)⋅∂𝜶⁡\(𝐱θ\)∂𝐱θ⋅∂𝐱θ∂θ,\\frac\{\\partial F\}\{\\partial\\theta\}=\\frac\{\\partial F\}\{\\partial\\bm\{\\alpha\}\(\\mathbf\{x\}\_\{\\theta\}\)\}\\cdot\\frac\{\\partial\\bm\{\\alpha\}\(\\mathbf\{x\}\_\{\\theta\}\)\}\{\\partial\\mathbf\{x\}\_\{\\theta\}\}\\cdot\\frac\{\\partial\\mathbf\{x\}\_\{\\theta\}\}\{\\partial\\theta\},where∂F∂αt​\(𝐱θ\)=f⁡\(𝐯t\)\\frac\{\\partial F\}\{\\partial\{\\alpha\}\_\{t\}\(\\mathbf\{x\}\_\{\\theta\}\)\}=f\(\\mathbf\{v\}\_\{t\}\)and∂𝜶⁡\(𝐱θ\)∂𝐱θ\\frac\{\\partial\\bm\{\\alpha\}\(\\mathbf\{x\}\_\{\\theta\}\)\}\{\\partial\\mathbf\{x\}\_\{\\theta\}\}comes from differentiating the steps our FW decomposition algorithm, Algorithm[1](https://arxiv.org/html/2609.25728#alg1)\. Note that under the stable\-oracle decision assumption introduced in Section[4\.1](https://arxiv.org/html/2609.25728#S4.SS1), the selected vertices𝐯t\\mathbf\{v\}\_\{t\}are locally constant functions of𝐱θ\\mathbf\{x\}\_\{\\theta\}almost everywhere\. Hence,∂𝐯t∂𝐱θ=0\\frac\{\\partial\\mathbf\{v\}\_\{t\}\}\{\\partial\\mathbf\{x\}\_\{\\theta\}\}=0wherever the derivative exists, and the gradient propagates only through the coefficients\. The gradient of the termλ​‖𝐱θ−𝐮⁡\(𝐱θ\)‖22\\lambda\\\|\\mathbf\{x\}\_\{\\theta\}\-\\mathbf\{u\}\(\\mathbf\{x\}\_\{\\theta\}\)\\\|\_\{2\}^\{2\}is computed similarly by differentiating through𝐮⁡\(𝐱θ\)=∑t=0T−1αt​\(𝐱θ\)​𝐯t\\mathbf\{u\}\(\\mathbf\{x\}\_\{\\theta\}\)=\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\(\\mathbf\{x\}\_\{\\theta\}\)\\mathbf\{v\}\_\{t\}\.

In Section[4\.1](https://arxiv.org/html/2609.25728#S4.SS1), we present our algorithm, which given𝐱\\mathbf\{x\}computes𝐮⁡\(𝐱\)∈𝒫\\mathbf\{u\}\(\\mathbf\{x\}\)\\in\\mathcal\{P\}satisfying the properties in \([2](https://arxiv.org/html/2609.25728#S4.E2)\)\. Moreover, we show that the weightsαt​\(𝐱\)\\alpha\_\{t\}\(\\mathbf\{x\}\)returned by our algorithm are a\.e\. differentiable w\.r\.t\.𝐱\\mathbf\{x\}and provide a convergence rate inTT, the number of iterations of the algorithm\. For notational simplicity, throughout Section[4\.1](https://arxiv.org/html/2609.25728#S4.SS1)we suppress the dependence on the target point𝐱\\mathbf\{x\}, writing𝐮\\mathbf\{u\}andαt\\alpha\_\{t\}in place of𝐮⁡\(𝐱\)\\mathbf\{u\}\(\\mathbf\{x\}\)andαt​\(𝐱\)\\alpha\_\{t\}\(\\mathbf\{x\}\)\.

### 4\.1From NN Output to Polytope: Frank\-Wolfe Decomposition

In this section, we describe our algorithm presented in[Algorithm1](https://arxiv.org/html/2609.25728#alg1)\. Given𝐱\\mathbf\{x\}, define the following quadratic objective:

J⁡\(𝐮\)=12​‖𝐱−𝐮‖22,𝐮∈𝒫\.J\(\\mathbf\{u\}\)=\\tfrac\{1\}\{2\}\\\|\\mathbf\{x\}\-\\mathbf\{u\}\\\|\_\{2\}^\{2\},\\qquad\\mathbf\{u\}\\in\\mathcal\{P\}\.\(5\)Starting at an initial point𝐮0\\mathbf\{u\}\_\{0\}in the polytope our goal is to iteratively move towards𝐮1,𝐮2,\.\.,𝐮T−1\\mathbf\{u\}\_\{1\},\\mathbf\{u\}\_\{2\},\.\.,\\mathbf\{u\}\_\{T\-1\}so that at each iteration we decreaseJ⁡\(𝐮t\)J\(\\mathbf\{u\}\_\{t\}\), and at the end𝐮T−1\\mathbf\{u\}\_\{T\-1\}has the desired aforementioned properties\. More specifically, we start with an initial point𝐮0\\mathbf\{u\}\_\{0\}which is a corner of the polytope𝒫\\mathcal\{P\}, i\.e\.𝐯0=𝐮0∈ext⁡\(𝒫\)\\mathbf\{v\}\_\{0\}=\\mathbf\{u\}\_\{0\}\\in\\mathrm\{ext\}\(\\mathcal\{P\}\)\. Then, in the first iteration, we find another corner of the polytope𝒫\\mathcal\{P\}, say𝐯1\\mathbf\{v\}\_\{1\}and set𝐮1=\(1−γ1\)​𝐮0\+γ1​𝐯1\\mathbf\{u\}\_\{1\}=\(1\-\\gamma\_\{1\}\)\\mathbf\{u\}\_\{0\}\+\\gamma\_\{1\}\\mathbf\{v\}\_\{1\}with an appropriate valueγ1∈\[0,1\]\\gamma\_\{1\}\\in\[0,1\], to be determined later\. By construction,𝐮1∈𝒫\\mathbf\{u\}\_\{1\}\\in\\mathcal\{P\}as it is a convex combination of two corners of the polytope\. This iterative process is repeated and at every iterationttwe have

𝐮t=\(1−γt\)​𝐮t−1\+γt​𝐯t\\displaystyle\\mathbf\{u\}\_\{t\}=\(1\-\\gamma\_\{t\}\)\\mathbf\{u\}\_\{t\-1\}\+\\gamma\_\{t\}\\mathbf\{v\}\_\{t\}\(6\)with an appropriate value forγt∈\[0,1\]\\gamma\_\{t\}\\in\[0,1\]and𝐯t∈ext⁡\(𝒫\)\\mathbf\{v\}\_\{t\}\\in\\mathrm\{ext\}\(\\mathcal\{P\}\)\. The algorithm terminates afterT−1T\-1iterations and returns𝐮T−1\\mathbf\{u\}\_\{T\-1\}, which will lie in𝒫\\mathcal\{P\}by construction\. At this point,𝐮T−1\\mathbf\{u\}\_\{T\-1\}is a convex combination of𝐮0=𝐯0,…,𝐯T−1\\mathbf\{u\}\_\{0\}=\\mathbf\{v\}\_\{0\},\\dots,\\mathbf\{v\}\_\{T\-1\}with weightsα0,…,αT−1\\alpha\_\{0\},\\dots,\\alpha\_\{T\-1\}that are calculated recursively with respect toγt\\gamma\_\{t\}s\. Specifically,∑t=0T−1αt=1\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}=1, withα0=∏i=1T−1\(1−γi\)\\alpha\_\{0\}=\\prod\_\{i=1\}^\{T\-1\}\(1\-\\gamma\_\{i\}\)andαt=γt​∏i=t\+1T−1\(1−γi\)\\alpha\_\{t\}=\\gamma\_\{t\}\\prod\_\{i=t\+1\}^\{T\-1\}\(1\-\\gamma\_\{i\}\)\.

Algorithm 1FW Decomposition\(𝐱,𝒫,T,LMO\)\(\\mathbf\{x\},\\mathcal\{P\},T;\\mathrm\{LMO\}\)1:Input:Target

𝐱∈ℝn\\mathbf\{x\}\\in\\mathbb\{R\}^\{n\}, polytope

𝒫\\mathcal\{P\}, support budget

TT, and an exact or approximate LMO\.

2:Initialize:

𝐯0=𝐮0=arg​max𝐯∈𝒫⁡⟨𝐱,𝐯⟩\\mathbf\{v\}\_\{0\}=\\mathbf\{u\}\_\{0\}=\\operatorname\*\{arg\\,max\}\_\{\\mathbf\{v\}\\in\\mathcal\{P\}\}\\langle\\mathbf\{x\},\\mathbf\{v\}\\rangle,

α0←1\\alpha\_\{0\}\\leftarrow 1\.

3:for

t=1t=1to

T−1T\-1do

4:

𝐫t←𝐱−𝐮t−1\\mathbf\{r\}\_\{t\}\\leftarrow\\mathbf\{x\}\-\\mathbf\{u\}\_\{t\-1\}\{Residual\}

5:

𝐯t←arg​max𝐯∈𝒫⁡⟨𝐫t,𝐯⟩\\mathbf\{v\}\_\{t\}\\leftarrow\\operatorname\*\{arg\\,max\}\_\{\\mathbf\{v\}\\in\\mathcal\{P\}\}\\langle\\mathbf\{r\}\_\{t\},\\mathbf\{v\}\\rangle\{Find a corner via LMO\}

6:

𝐝t←𝐯t−𝐮t−1\\mathbf\{d\}\_\{t\}\\leftarrow\\mathbf\{v\}\_\{t\}\-\\mathbf\{u\}\_\{t\-1\}\{Descent direction\}

7:

γt←Clip⁡\(0,1,⟨𝐫t,𝐝t⟩/‖𝐝t‖22\)\\gamma\_\{t\}\\leftarrow\\mathrm\{Clip\}\(0,1,\\langle\\mathbf\{r\}\_\{t\},\\mathbf\{d\}\_\{t\}\\rangle/\\\|\\mathbf\{d\}\_\{t\}\\\|\_\{2\}^\{2\}\)\{Line search step\}

8:

αj←\(1−γt\)​αj\\alpha\_\{j\}\\leftarrow\(1\-\\gamma\_\{t\}\)\\alpha\_\{j\}for all existing weight

αj\\alpha\_\{j\},

j∈\{0,…,t−1\}j\\in\\\{0,\\dots,t\-1\\\}\.

9:

αt←γt\\alpha\_\{t\}\\leftarrow\\gamma\_\{t\}
10:

𝐮t←\(1−γt\)​𝐮t−1\+γt​𝐯t\\mathbf\{u\}\_\{t\}\\leftarrow\(1\-\\gamma\_\{t\}\)\\mathbf\{u\}\_\{t\-1\}\+\\gamma\_\{t\}\\mathbf\{v\}\_\{t\}
11:endfor

12:Output:Weights and vertices

\{\(αt,𝐯t\)\}t=0T−1\\\{\(\\alpha\_\{t\},\\mathbf\{v\}\_\{t\}\)\\\}\_\{t=0\}^\{T\-1\}such that

𝐮T−1=∑t=0T−1αt​𝐯t\\mathbf\{u\}\_\{T\-1\}=\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\\mathbf\{v\}\_\{t\}\.

The main two points to address here are which corner𝐯t\\mathbf\{v\}\_\{t\}to choose and what is a good value for step sizeγt\\gamma\_\{t\}\.

Note that the functionJ⁡\(𝐮\)J\(\\mathbf\{u\}\)is convex and its gradient is∇J​\(𝐮\)=𝐮−𝐱\\nabla J\(\\mathbf\{u\}\)=\\mathbf\{u\}\-\\mathbf\{x\}\. We define theresidualat iterationttas𝐫t:=𝐱−𝐮t−1=−∇J​\(𝐮t−1\)\\mathbf\{r\}\_\{t\}:=\\mathbf\{x\}\-\\mathbf\{u\}\_\{t\-1\}=\-\\nabla J\(\\mathbf\{u\}\_\{t\-1\}\), which points from our current estimate toward the target\. Frank–Wolfe is a "projection\-free" method\. At every iteration, instead of moving toward the target in the gradient direction and then projecting back into𝒫\\mathcal\{P\}\(which can be computationally expensive\), it finds the vertex𝐯t\\mathbf\{v\}\_\{t\}of the polytope that lies furthest in the direction of the residual𝐫t\\mathbf\{r\}\_\{t\}\. It then takes a small stepγt\\gamma\_\{t\}from the current point𝐮t−1\\mathbf\{u\}\_\{t\-1\}toward that vertex𝐯t\\mathbf\{v\}\_\{t\}\. Since both𝐮t−1\\mathbf\{u\}\_\{t\-1\}and𝐯t\\mathbf\{v\}\_\{t\}are in𝒫\\mathcal\{P\}, their convex combination is guaranteed to remain inside𝒫\\mathcal\{P\}\. The exact procedure is shown in[Algorithm1](https://arxiv.org/html/2609.25728#alg1), and in order to find𝐯t\\mathbf\{v\}\_\{t\}we require what is known as LMO\.

###### Definition 4\.2\(Linear Maximization Oracle \(LMO\)\)\.

Given a residual𝐫∈ℝn\\mathbf\{r\}\\in\\mathbb\{R\}^\{n\}, an exact LMO returns a vertex𝐯⋆\\mathbf\{v\}^\{\\star\}that maximizes the alignment with𝐫\\mathbf\{r\}:

𝐯⋆=arg​max⁡⟨𝐫,𝐯⟩s\.t\.𝐯∈𝒫\\displaystyle\\mathbf\{v\}^\{\\star\}=\\operatorname\*\{arg\\,max\}\\langle\\mathbf\{r\},\\mathbf\{v\}\\rangle\\quad\\text\{s\.t\.\}\\quad\\mathbf\{v\}\\in\\mathcal\{P\}\(7\)

A fundamental property of linear optimization is that a linear objective over a polytope attains an optimum at a vertex\. Hence, an LMO returns a \(integral\) vertex of the polytope\. When exact maximization is difficult, the convergence statement below uses a relative Frank–Wolfe gap approximation\[[21](https://arxiv.org/html/2609.25728#bib.bib3),[34](https://arxiv.org/html/2609.25728#bib.bib8),[42](https://arxiv.org/html/2609.25728#bib.bib51),[70](https://arxiv.org/html/2609.25728#bib.bib4)\]\. Formally, at iterationtt, aδ\\delta\-relative FW\-gap oracle returns𝐯t∈ext⁡\(𝒫\)\\mathbf\{v\}\_\{t\}\\in\\mathrm\{ext\}\(\\mathcal\{P\}\)satisfying

⟨𝐫t,𝐯t−𝐮t−1⟩≥δ​max𝐯∈𝒫​⟨𝐫t,𝐯−𝐮t−1⟩,\\displaystyle\\langle\\mathbf\{r\}\_\{t\},\\mathbf\{v\}\_\{t\}\-\\mathbf\{u\}\_\{t\-1\}\\rangle\\;\\geq\\;\\delta\\max\_\{\\mathbf\{v\}\\in\\mathcal\{P\}\}\\langle\\mathbf\{r\}\_\{t\},\\mathbf\{v\}\-\\mathbf\{u\}\_\{t\-1\}\\rangle,\(8\)where𝐫t=𝐱−𝐮t−1\\mathbf\{r\}\_\{t\}=\\mathbf\{x\}\-\\mathbf\{u\}\_\{t\-1\}andδ∈\(0,1\]\\delta\\in\(0,1\]\. This is the approximation notion used in our convergence statement\. Exact LMOs satisfy \([8](https://arxiv.org/html/2609.25728#S4.E8)\) withδ=1\\delta=1\.

For the step sizeγt\\gamma\_\{t\}, several choices are possible\. We compute it using a simple line\-search:

γt=Clip⁡\(0,1,⟨𝐫t,𝐝t⟩/‖𝐝t‖22\),\\gamma\_\{t\}=\\mathrm\{Clip\}\(0,1,\\langle\\mathbf\{r\}\_\{t\},\\mathbf\{d\}\_\{t\}\\rangle/\\\|\\mathbf\{d\}\_\{t\}\\\|\_\{2\}^\{2\}\),with the conventionγt=0\\gamma\_\{t\}=0when𝐝t=0\\mathbf\{d\}\_\{t\}=0\. This allows for efficient computation and supports branchwise differentiation of the weights\{αt\}\\\{\\alpha\_\{t\}\\\}\.

Assumption \(stable oracle decisions\)\.Fix the support budgetTT\. We assume that the deterministic oracle and tie\-breaking rule do not switch pathologically\. The set of inputs𝐱\\mathbf\{x\}at which an arbitrarily small perturbation can change one of the selected vertices, or change whether a line\-search step is clipped, has zero volume\. Equivalently, for almost every input𝐱\\mathbf\{x\}, there is a small neighborhood around𝐱\\mathbf\{x\}on which Algorithm[1](https://arxiv.org/html/2609.25728#alg1)selects the same vertices𝐯0,…,𝐯T−1\\mathbf\{v\}\_\{0\},\\ldots,\\mathbf\{v\}\_\{T\-1\}and uses the same line\-search clipping cases\. This assumption rules out artificial rules that switch between equally valid approximate vertices on dense sets\. It is satisfied by the standard comparison\-based oracles used in combinatorial optimization, including top\-kkselection, matroid greedy algorithms such as Kruskal’s algorithm, and exact assignment or matching oracles with fixed tie\-breaking\.

All these results together yield the following theorem\. The convergence guarantee is the standard exact Frank–Wolfe rate, and the a\.e\. differentiability proof is given in[AppendixA](https://arxiv.org/html/2609.25728#A1)\.

###### Theorem 4\.3\(Exact Frank–Wolfe Decomposition\)\.

Let𝐱∈ℝn\\mathbf\{x\}\\in\\mathbb\{R\}^\{n\}, and let𝒫⊂ℝn\\mathcal\{P\}\\subset\\mathbb\{R\}^\{n\}be a polytope with diameterD:=diam⁡\(𝒫\)=max𝐮,𝐯∈𝒫⁡‖𝐮−𝐯‖2D:=\\mathrm\{diam\}\(\\mathcal\{P\}\)=\\max\_\{\\mathbf\{u\},\\mathbf\{v\}\\in\\mathcal\{P\}\}\\\|\\mathbf\{u\}\-\\mathbf\{v\}\\\|\_\{2\}\. Run Algorithm[1](https://arxiv.org/html/2609.25728#alg1)with an exact LMO and support budgetT≥2T\\geq 2, so that it performsT−1T\-1Frank–Wolfe updates\. Then the algorithm outputs𝐮T−1=∑t=0T−1αt​𝐯t∈𝒫\\mathbf\{u\}\_\{T\-1\}=\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\\mathbf\{v\}\_\{t\}\\in\\mathcal\{P\}, where𝐯t∈ext⁡\(𝒫\)\\mathbf\{v\}\_\{t\}\\in\\operatorname\{ext\}\(\\mathcal\{P\}\),αt≥0\\alpha\_\{t\}\\geq 0, and∑t=0T−1αt=1\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}=1\. Moreover, if𝐮⋆=arg​min𝐮∈𝒫⁡J​\(𝐮\)\\mathbf\{u\}^\{\\star\}=\\operatorname\*\{arg\\,min\}\_\{\\mathbf\{u\}\\in\\mathcal\{P\}\}J\(\\mathbf\{u\}\), then

J⁡\(𝐮T−1\)−J⁡\(𝐮⋆\)≤2​D2T\+1\.J\(\\mathbf\{u\}\_\{T\-1\}\)\-J\(\\mathbf\{u\}^\{\\star\}\)\\;\\leq\\;\\frac\{2D^\{2\}\}\{T\+1\}\.The coefficientsαt\\alpha\_\{t\}are almost everywhere differentiable functions of𝐱\\mathbf\{x\}\.

###### Corollary 4\.4\(Approximate Frank–Wolfe Decomposition\)\.

If the exact LMO in[Theorem4\.3](https://arxiv.org/html/2609.25728#S4.Thmtheorem3)is replaced by aδ\\delta\-relative FW\-gap oracle satisfying \([8](https://arxiv.org/html/2609.25728#S4.E8)\), then the same convex\-combination conclusion holds\. LetD=diam⁡\(𝒫\)D=\\mathrm\{diam\}\(\\mathcal\{P\}\),𝐮⋆=arg​min𝐮∈𝒫⁡J​\(𝐮\)\\mathbf\{u\}^\{\\star\}=\\operatorname\*\{arg\\,min\}\_\{\\mathbf\{u\}\\in\\mathcal\{P\}\}J\(\\mathbf\{u\}\), andh0=J⁡\(𝐮0\)−J⁡\(𝐮⋆\)h\_\{0\}=J\(\\mathbf\{u\}\_\{0\}\)\-J\(\\mathbf\{u\}^\{\\star\}\)\. Then

J⁡\(𝐮T−1\)−J⁡\(𝐮⋆\)≤2​\(D2/δ\+h0\)δ⁡\(T−1\)\+2\.J\(\\mathbf\{u\}\_\{T\-1\}\)\-J\(\\mathbf\{u\}^\{\\star\}\)\\;\\leq\\;\\frac\{2\\left\(D^\{2\}/\\delta\+h\_\{0\}\\right\)\}\{\\delta\(T\-1\)\+2\}\.The coefficient map remains almost everywhere differentiable\.

Discussion on projection onto the polytope\.In the above discussion we did not require neural network output to lie in the polytope, yielding a generic projection\-free framework\. Nevertheless, when this occurs—either by design of the neural network or because it is empirically advantageous—our framework aligns with*approximate Carathéodory*results\. When𝐱∈𝒫\\mathbf\{x\}\\in\\mathcal\{P\}, minimizingJ⁡\(𝐮\)=12​‖𝐱−𝐮‖22J\(\\mathbf\{u\}\)=\\tfrac\{1\}\{2\}\\\|\\mathbf\{x\}\-\\mathbf\{u\}\\\|\_\{2\}^\{2\}via[Algorithm1](https://arxiv.org/html/2609.25728#alg1)produces a sparse convex combination of vertices of𝒫\\mathcal\{P\}approximating𝐱\\mathbf\{x\}\. Constructing such sparse convex combinations is the goal of approximate Carathéodory results\. Unlike the exact Carathéodory theorem, the approximate variants relax exact representation in favor of an additive error to achieve sparser convex combinations\[[51](https://arxiv.org/html/2609.25728#bib.bib23),[4](https://arxiv.org/html/2609.25728#bib.bib12),[46](https://arxiv.org/html/2609.25728#bib.bib1)\]\. In the exact\-LMO setting,[Theorem4\.3](https://arxiv.org/html/2609.25728#S4.Thmtheorem3)yields the approximate Carathéodory bound‖𝐱−𝐮T−1‖22≤4​D2T\+1\.\\\|\\mathbf\{x\}\-\\mathbf\{u\}\_\{T\-1\}\\\|\_\{2\}^\{2\}\\leq\\tfrac\{4D^\{2\}\}\{T\+1\}\.

### 4\.2Examples with Efficient \(Approximate\) LMO

We briefly discuss how our framework extends to structured polytopes; details are in[AppendixG](https://arxiv.org/html/2609.25728#A7)\.

Matroid polytopes\.Matroid base polytopes capture many combinatorial constraints \(e\.g\., cardinality, partition, spanning trees\) and admit efficient LMOs via greedy maximum\-weight base algorithms\. Since our method requires only an efficient LMO, it applies uniformly across matroids without modification, unlike\[[23](https://arxiv.org/html/2609.25728#bib.bib11)\], which requires constraint\-specific projections and decompositions\.

Birkhoff and matching polytopes\.The Birkhoff and \(perfect\) matching polytopes admit LMOs via maximum\-weight matching, allowing direct optimization without projections or penalties\. This is unlike prior work \(e\.g\.,\[[50](https://arxiv.org/html/2609.25728#bib.bib21)\]\), which relies on specialized decompositions and projection\. Although exact LMOs may be costly at scale, we find that simple greedy approximations suffice empirically\.

## 5Applications and Experiments

We evaluate our framework on three core CO problems—MC, QAP, and TSP—and show that the same decomposition and extension apply across all cases\. Our method is competitive with and often outperforms SOTA approaches\. Compared to SOTA approaches, which rely on problem\-specific and novel engineered constructions, this generality is a key advantage\. For testing methods and choice of baselines, we keep a strict one\-shot setting \(e\.g\., without Test Time Optimization \(TTO\)\) to ensure fairness, detailed in[AppendixE](https://arxiv.org/html/2609.25728#A5)\. Code and data are available in[link](https://anonymous.4open.science/r/FW_NCO-FA5C)\.

### 5\.1Maximum Coverage

Table 1:MC performance comparison across three datasets, reporting mean inference time and coverage acrosskk\. Among neural methods,■\\blacksquareindicates ranking the 1st,■\\blacksquarethe 2nd in each column\.Let𝒰\\mathcal\{U\}be a ground set of elements with weightswu≥0w\_\{u\}\\geq 0foru∈𝒰u\\in\\mathcal\{U\}, and let\{S1,…,Sn\}\\\{S\_\{1\},\\dots,S\_\{n\}\\\}be subsetsSi⊆𝒰S\_\{i\}\\subseteq\\mathcal\{U\}\. The value of𝐱∈\{0,1\}n\\mathbf\{x\}\\in\\\{0,1\\\}^\{n\}is the total weight of covered elements,fcov\(𝐱\)=∑u∈𝒰wu\[∃i∈\[n\]withu∈Siand𝐱\(i\)=1\]\.f\_\{\\text\{cov\}\}\(\\mathbf\{x\}\)\\;=\\;\\sum\_\{u\\in\\mathcal\{U\}\}w\_\{u\}\\,\\mathbf\{1\}\\\!\\left\[\\exists\\,i\\in\[n\]\\ \\text\{with\}\\ u\\in S\_\{i\}\\ \\text\{and\}\\ \\mathbf\{x\}\(i\)=1\\right\]\.The MC problem under cardinality constraint ismax𝐱∈\{0,1\}n⁡fcov​\(𝐱\)\\max\_\{\\mathbf\{x\}\\in\\\{0,1\\\}^\{n\}\}f\_\{\\text\{cov\}\}\(\\mathbf\{x\}\)subject to∑i=1n𝐱⁡\(i\)=k\\sum\_\{i=1\}^\{n\}\\mathbf\{x\}\(i\)=k, wherekkdenotes the prescribed cardinality\.

Neural model and continuous targets\.This problem can naturally be represented as a bipartite graph where the goal is to selectkknodes from one part so that we maximize the total number of their neighbors in the other part\. We apply a GAT\-based encoder to obtain node embeddings\. This𝐱θ\\mathbf\{x\}\_\{\\theta\}is the input to ourFW Decompositionalgorithm to obtain a small list of subsets\{S\(t\)\}t=0T−1\\\{S^\{\(t\)\}\\\}\_\{t=0\}^\{T\-1\}, of sizekk, and coefficients\{αt\}t=0T−1\\\{\\alpha\_\{t\}\\\}\_\{t=0\}^\{T\-1\}\. For training, we minimize the differentiable surrogate lossℒCOV=−∑t=0T−1αtfcov\(𝟏S\(t\)\)\+λ∥𝐱θ−∑t=0T−1αt𝟏S\(t\)∥22\\mathcal\{L\}\_\{\\text\{COV\}\}=\-\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\\,f\_\{\\text\{cov\}\}\(\\mathbf\{1\}\_\{S^\{\(t\)\}\}\)\+\\lambda\\\|\\mathbf\{x\}\_\{\\theta\}\-\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\\mathbf\{1\}\_\{S^\{\(t\)\}\}\\\|\_\{2\}^\{2\}and backpropagate through the FW steps into the GAT parameters\. Following\[[23](https://arxiv.org/html/2609.25728#bib.bib11)\]\(GeoNCO\), we train on synthetic random graphs and test on real\-world graphs\. For evaluation,[Table1](https://arxiv.org/html/2609.25728#S5.T1), we reportmaxt⁡fcov​\(𝟏S\(t\)\)\\max\_\{t\}f\_\{\\text\{cov\}\}\(\\mathbf\{1\}\_\{S^\{\(t\)\}\}\)without TTO and compare with the baselines without TTO\. Our method is faster thanGeoNCOand it often outperforms the learning baselines\. Specifically, it consistently stays on the pareto front over the learning baselines and is also competitive against non\-learning baselines considering the quality\-efficiency tradeoff\. \(See Section[B\.4](https://arxiv.org/html/2609.25728#A2.SS4)for adversarial instances where greedy fails\.\)

Ablation on projection for MC\.[Table4](https://arxiv.org/html/2609.25728#A2.T4)shows that removing the polytope projection preserves solution quality while improving runtime\. Although the model converges to a point outside the hypersimplex, the solution quality from our decomposition is slightly better than the projected version\. This suggests the projection\-free variant allows greater exploration by avoiding hard constraints\.

### 5\.2Quadratic Assignment Problem

There arennfacilities andnnlocations\. Distances between locations are given by a matrix𝐁∈ℝn×n\\mathbf\{B\}\\in\\mathbb\{R\}^\{n\\times n\}, and flows between facilities by a matrix𝐀∈ℝn×n\\mathbf\{A\}\\in\\mathbb\{R\}^\{n\\times n\}\. The goal is to assign each facility to a distinct location so as to minimize the total flow\-weighted distance\. LetΠn\\Pi\_\{n\}denote the set ofn×nn\\times npermutation matrices\. The objective ismin𝐏∈Πn⁡fqap​\(𝐏\)\\min\_\{\\mathbf\{P\}\\in\\Pi\_\{n\}\}f\_\{\\text\{qap\}\}\(\\mathbf\{P\}\)wherefqap​\(𝐏\)=trace⁡\(𝐀𝐏𝐁𝐏⊤\)f\_\{\\text\{qap\}\}\(\\mathbf\{P\}\)=\\mathrm\{trace\}\(\\mathbf\{A\}\\mathbf\{P\}\\mathbf\{B\}\\mathbf\{P\}^\{\\top\}\)for all𝐏∈Πn\\mathbf\{P\}\\in\\Pi\_\{n\}\.

Neural model and continuous targets\.Given the output of encoder network, GraphSAGE,𝐗θ∈\[0,1\]n×n\\mathbf\{X\}\_\{\\theta\}\\in\[0,1\]^\{n\\times n\}, we apply ourFW Decompositionalgorithm using a deterministic greedy algorithm for weighted maximum matching\. AfterTTiterations, it returns a short list of permutation matrices\{𝐏t\}t=0T−1\\\{\\mathbf\{P\}\_\{t\}\\\}\_\{t=0\}^\{T\-1\}and coefficients\{αt\}t=0T−1\\\{\\alpha\_\{t\}\\\}\_\{t=0\}^\{T\-1\}\. We train directly against the discrete QAP objective and using the followingℒQAP​\(𝐗θ\)=∑t=0T−1αt​fqap​\(𝐏t\)\+λ​‖𝐗θ−∑t=0T−1αt​𝐏t‖F2\\mathcal\{L\}\_\{\\text\{QAP\}\}\(\\mathbf\{X\}\_\{\\theta\}\)=\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\\,f\_\{\\text\{qap\}\}\(\\mathbf\{P\}\_\{t\}\)\+\\lambda\\\|\\mathbf\{X\}\_\{\\theta\}\-\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\\,\\mathbf\{P\}\_\{t\}\\\|\_\{F\}^\{2\}\. Gradients are backpropagated through the FW steps into the parameters that produce𝐗θ\\mathbf\{X\}\_\{\\theta\}\. At evaluation we reportmint⁡fqap​\(𝐏t\)\\min\_\{t\}f\_\{\\text\{qap\}\}\(\\mathbf\{P\}\_\{t\}\)\.

We train on synthetic instances following\[[59](https://arxiv.org/html/2609.25728#bib.bib52)\]\. For sizenn, we sample locations𝐋∼Uniform​\(0,1\)2\\mathbf\{L\}\\sim\\mathrm\{Uniform\}\(0,1\)^\{2\}and a symmetric flow matrix𝐀∈ℝn×n\\mathbf\{A\}\\in\\mathbb\{R\}^\{n\\times n\}with zero diagonal, upper\-triangular entries sampled from\[0,1\]\[0,1\], and entries independently set to zero with probabilitypp\(symmetrically\)\. Testing uses the standardQAPLIBbenchmark\[[8](https://arxiv.org/html/2609.25728#bib.bib61)\]\(134 instances, 15 classes; see[SectionC\.1](https://arxiv.org/html/2609.25728#A3.SS1)\)\.

Consistent with prior studies, we report per\-class and aggregated gap statistics on instances of size 12–64, where baseline methods are tractable\. Recent learning\-based QAP methods, such asRGM\[[41](https://arxiv.org/html/2609.25728#bib.bib53)\]andSAWT\[[59](https://arxiv.org/html/2609.25728#bib.bib52)\], rely on iterative Learn\-to\-Construct or Learn\-to\-Improve pipelines which could be time consuming and prevent scalability\. Our approach, on the other hand, produces a permutation in a single forward pass, which compared to the learning baselines, yields better efficiency, while achieving a better out of distribution generalization \(OOD\) quality\.

Table 2:Mean gaps and inference time \(s\) for models pre\-trained on QAP32 and applied to QAPLIB instances\.Std\.is the standard deviation of the class\-wise mean gaps across the 14 data classes \(bur–wil\)\. In each column,■\\blacksquareindicates ranking the 1st,■\\blacksquarethe 2nd\.In our experiments,[Table2](https://arxiv.org/html/2609.25728#S5.T2), our method achieves the best overall performance among all baselines, including both learning\-based and classical solvers\. The consistently low standard deviation of the average optimality gap across instance classes indicates stable behavior and good generalization across distributions\. Moreover, our approach is the fastest neural method, running noticeably faster than other neural baselines\. The only two methods with lower runtime,SMandRRWM, are purely traditional heuristics, which come at a substantial loss in solution quality\. While[Table2](https://arxiv.org/html/2609.25728#S5.T2)reports sizes 12–64, our method scales to larger instances \(\>64\>64\) with similar average gaps \([Table8](https://arxiv.org/html/2609.25728#A3.T8)\)\. Most prior methods do not scale to such sizes with reasonable compute resources\.

Ablation on projection for QAP\.[Algorithm1](https://arxiv.org/html/2609.25728#alg1)does not theoretically require its input to be in the Birkhoff polytope\. However, we observed a weaker OOD generalization results for QAP due to the strong heterogeneity and distribution shift present inQAPLIB\. Unconstrained outputs, i\.e\., "far from the polytope", cause the decomposition algorithm to return low quality permutations under drastic distributional shift\. We therefore include a lightweight Sinkhorn projection layer that enforces approximate doubly\-stochasticity, improving generalization\. This step is introduced for practical considerations, not theoretical necessity\.

### 5\.3TSP via Spanning Trees and Learned Matchings

We consider spanning tree polytopes, where𝒞\\mathcal\{C\}is the set of spanning trees of a graph\. For eachT∈𝒞T\\in\\mathcal\{C\}, the objectivefTSP​\(T\)f\_\{\\text\{TSP\}\}\(T\)is the cost of the TSP tour obtained deterministically via Christofides’s algorithm\[[11](https://arxiv.org/html/2609.25728#bib.bib70)\]\. Training learns a distribution over spanning trees optimized for this objective\.

Algorithm\-guided neural optimization for TSP\.One of the appeals of our approach is that it allows us to infuse learnable NN modules into an algorithmic framework\. We adopt an algorithmically aligned, data\-driven approach that integrates learning into the classical Christofides–Serdyukov algorithmic framework\. Christofides’s algorithm builds a TSP solution from a minimum spanning tree, whereas recent advances \(e\.g\.,\[[3](https://arxiv.org/html/2609.25728#bib.bib54),[25](https://arxiv.org/html/2609.25728#bib.bib56),[26](https://arxiv.org/html/2609.25728#bib.bib55)\]\) sample from spanning tree distributions constructed from a subtour LP solution before applying variants of the Christofides pipeline\. Our approach can be viewed as a learning\-based analogue of this paradigm: instead of constructing such distributions analytically, we learn a distribution over spanning trees in an end\-to\-end manner, optimized directly for the final TSP objective\. Moreover, we extend this idea further by incorporating a learning module for the matching step, effectively learning over matching polytopes as well\. Given the neural network output,FW Decompositionreturns a sparse distribution over spanning trees using Kruskal’s algorithm as LMO\. We sample trees from this distribution, and recover tours by applying minimum\-weight perfect matching on odd degree nodes followed by shortcutting\. This pipeline defines a differentiable surrogate objective that mirrors the Christofides algorithm end to end\. This also allows the perfect\-matching step itself to be learned jointly, improving solution quality\.

In our experiments, we consider Euclidean TSP instances\(G=\(V,E\),D\)\(G=\(V,E\),D\), whereVVis a set of points drawn uniformly at random from the unit square\[0,1\]2\[0,1\]^\{2\},EEis the complete edge set onVV, andDDis the corresponding symmetric distance matrix induced by Euclidean distances, which satisfies the triangle inequality\. A tour is defined as a Hamiltonian cycle that visits each node exactly once, and its cost is given by the sum of the lengths of its constituent edges underDD\. The objective is to find a tour of minimum total cost\.

Neural TSP methods commonly fall into two paradigms\[[44](https://arxiv.org/html/2609.25728#bib.bib19)\]: Learning\-to\-Construct \(LC\), which builds TSP tours sequentially with hard feasibility checks but suffers from high inference latency on large instances due to its inherently sequential decoding, and Global Prediction \(GP\), which predicts a full structure \(e\.g\., a heatmap\) and relies on post\-processing, yielding fast inference but strong sensitivity to recovery heuristics\. COExpander\[[44](https://arxiv.org/html/2609.25728#bib.bib19)\]combines LC and GP via Adaptive Expansion \(AE\), but still depends on heatmap recovery, which still is a combinatorial optimization subroutine, inheriting GP’s core limitation\. Motivated by critiques of heatmap\-guided post\-hoc search\[[68](https://arxiv.org/html/2609.25728#bib.bib64)\], our method directly optimizes feasible tours end\-to\-end by differentiating through the selected FW branch, reducing reliance on hand\-crafted post\-processing\.

We compare against four baseline categories: \(i\) exact solvers \(Concorde,Gurobi\); \(ii\) classical heuristics \(LKH3\[[20](https://arxiv.org/html/2609.25728#bib.bib66)\], 500 trials\); \(iii\) neural LC methods \(Sym\-NCO\[[29](https://arxiv.org/html/2609.25728#bib.bib69)\]\), implemented in the RL4CO library\[[5](https://arxiv.org/html/2609.25728#bib.bib67)\]; and \(iv\) neural GP/AE and heatmap\-based methods across RL and supervised settings, includingDIMES\[[52](https://arxiv.org/html/2609.25728#bib.bib58)\],DIFUSCO\[[58](https://arxiv.org/html/2609.25728#bib.bib57)\], andCOExpander, a SOTA AE method\. In[AppendixE](https://arxiv.org/html/2609.25728#A5), we provide discussion regarding approaches that use inference optimization techniques such as Gradient\-based Active Search, and Search\-based Refinement such as MCTS and 2\-opt\.

Table 3:Performance comparison on TSP benchmarks across different scales\. Results include optimality gap \(%\) and inference time\. Methods are categorized as Unsupervised Learning \(UL\), Exact, Heuristics, Reinforcement Learning \(RL\), and Supervised Learning \(SL\)\. Among neural methods,■\\blacksquareindicates ranking the 1st,■\\blacksquarethe 2nd in each column\.Our results,[Table3](https://arxiv.org/html/2609.25728#S5.T3), show that our method consistently outperforms unsupervised neural baselines including RL\. On larger instances, it also surpassesDIFUSCO, a strong diffusion\-based approach trained with supervision\. While we are still behind the strongest supervised baselines, e\.g\.COExpander, our approach is more resource\-efficient\. In particular,DIFUSCOandCOExpanderrely on up to∼\\sim1\.5M labeled instances \(depending on the instance size\) and require extensive training time on good hardware, whereas our model can be trained in under one hour on a single outdated GPU\.

Ablation for TSP\. We examine the effect of projecting the neural network output onto the spanning tree polytope and find that removing projection slightly improves performance on both TSP\-500 and TSP\-1000 \([Table13](https://arxiv.org/html/2609.25728#A4.T13)\)\. We also study generalization across instance sizes and observe strong cross\-size generalization \([Table12](https://arxiv.org/html/2609.25728#A4.T12)\)\. Our model uses subtour LP solutions\[[14](https://arxiv.org/html/2609.25728#bib.bib31)\]as features, obtained via the QSopt LP solver, also used by Concorde\. Despite this shared solver, our learning pipeline is significantly faster on large instances\. Our ablation study,[SectionD\.7](https://arxiv.org/html/2609.25728#A4.SS7), confirm that these gains arise from learned structures rather than arbitrary perturbations of the LP relaxation\.

## 6Conclusion

We propose a general framework for self\-supervised CO under constraints\. Our approach is flexible in that it does not require the neural network output to lie in a feasible polytope\. It uses a generic procedure, to decompose the output and obtain a sparse distribution over feasible solutions, effectively enabling SSL\. Our method shows strong empirical performance compared to neural baselines on several CO problems\. We present an example of infusing learned modules within an algorithm in an end\-to\-end pipeline, showcasing this approach on TSP using the Christofides algorithm\. We believe this paves the way toward new NN\-infused algorithms that can effectively tackle complex CO problems\. We note that model architecture can affect performance; we do not study this here and leave it for future work\.

## References

- \[1\]S\. Amizadeh, S\. Matusevych, and M\. Weimer\(2018\)Learning to solve circuit\-sat: an unsupervised differentiable approach\.InInternational Conference on Learning Representations \(ICLR\),Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[2\]D\. L\. Applegate, R\. E\. Bixby, V\. Chvátal, and W\. J\. Cook\(2011\)The traveling salesman problem: a computational study\.InThe Traveling Salesman Problem,Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p1.1)\.
- \[3\]A\. Asadpour, M\. X\. Goemans, A\. Madry, S\. O\. Gharan, and A\. Saberi\(2017\)An o \(log n/log log n\)\-approximation algorithm for the asymmetric traveling salesman problem\.Operations Research65\(4\),pp\. 1043–1061\.Cited by:[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p2.1)\.
- \[4\]S\. Barman\(2015\)Approximating nash equilibria and dense bipartite subgraphs via an approximate version of caratheodory’s theorem\.InProceedings of the forty\-seventh annual ACM symposium on Theory of computing,pp\. 361–369\.Cited by:[§1](https://arxiv.org/html/2609.25728#S1.p4.1),[§2](https://arxiv.org/html/2609.25728#S2.p3.1),[§4\.1](https://arxiv.org/html/2609.25728#S4.SS1.p8.1)\.
- \[5\]F\. Berto, C\. Hua, J\. Park, L\. Luttmann, Y\. Ma, F\. Bu, J\. Wang, H\. Ye, M\. Kim, S\. Choi,et al\.\(2025\)Rl4co: an extensive reinforcement learning for combinatorial optimization benchmark\.InProceedings of the 31st ACM SIGKDD Conference on Knowledge Discovery and Data Mining V\. 2,pp\. 5278–5289\.Cited by:[§D\.3](https://arxiv.org/html/2609.25728#A4.SS3.p4.1),[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p5.1)\.
- \[6\]M\. Böther, O\. Kißig, M\. Taraz, S\. Cohen, K\. Seidel, and T\. Friedrich\(2022\)What’s wrong with deep learning in tree search for combinatorial optimization\.InInternational Conference on Learning Representations \(ICLR\),Cited by:[Appendix E](https://arxiv.org/html/2609.25728#A5.SS0.SSS0.Px1.p2.1)\.
- \[7\]F\. Bu, H\. Jo, S\. Y\. Lee, S\. Ahn, and K\. Shin\(2024\)Tackling prevalent conditions in unsupervised combinatorial optimization: cardinality, minimum, covering, and more\.In41st International Conference on Machine Learning \(ICML\),pp\. 4696–4729\.Cited by:[§B\.1](https://arxiv.org/html/2609.25728#A2.SS1.SSS0.Px1.p1.1),[§B\.2](https://arxiv.org/html/2609.25728#A2.SS2.p3.1),[§B\.2](https://arxiv.org/html/2609.25728#A2.SS2.p4.1),[Appendix E](https://arxiv.org/html/2609.25728#A5.SS0.SSS0.Px1.p2.1),[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[8\]R\. E\. Burkard, S\. E\. Karisch, and F\. Rendl\(1997\)QAPLIB–a quadratic assignment problem library\.Journal of Global optimization10\(4\),pp\. 391–403\.Cited by:[§C\.1](https://arxiv.org/html/2609.25728#A3.SS1.SSS0.Px1.p2.1),[§5\.2](https://arxiv.org/html/2609.25728#S5.SS2.p3.1)\.
- \[9\]M\. D\. Canon and C\. D\. Cullum\(1968\)A tight upper bound on the rate of convergence of Frank–Wolfe algorithm\.SIAM Journal on Control6\(4\),pp\. 509–516\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p3.1)\.
- \[10\]C\. Carathéodory\(1907\)Über den variabilitätsbereich der koeffizienten von potenzreihen, die gegebene werte nicht annehmen\.Mathematische Annalen64\(1\),pp\. 95–115\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p3.1)\.
- \[11\]N\. Christofides\(1976\)Worst\-case analysis of a new heuristic for the travelling salesman problem\.Technical reportCited by:[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p1.1)\.
- \[12\]C\. W\. Combettes and S\. Pokutta\(2023\)Revisiting the approximate Carathéodory problem via the Frank\-Wolfe algorithm\.Mathematical Programming197\(1\),pp\. 191–214\.External Links:ISSN 1436\-4646Cited by:[§1](https://arxiv.org/html/2609.25728#S1.p4.1),[§2](https://arxiv.org/html/2609.25728#S2.p3.1)\.
- \[13\]Y\. Crama\(1993\)Concave extensions for nonlinear 0–1 maximization problems\.Mathematical Programming61\(1\),pp\. 53–60\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p1.1)\.
- \[14\]G\. Dantzig, R\. Fulkerson, and S\. Johnson\(1954\)Solution of a large\-scale traveling\-salesman problem\.Journal of the operations research society of America2\(4\),pp\. 393–410\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p1.1),[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p7.1)\.
- \[15\]P\. L\. Donti, D\. Rolnick, and J\. Z\. Kolter\(2021\)DC3: a learning method for optimization with hard constraints\.InInternational Conference on Learning Representations \(ICLR\),Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[16\]J\. E\. Falk and K\. R\. Hoffman\(1976\)A successive underestimation method for concave minimization problems\.Mathematics of operations research1\(3\),pp\. 251–259\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p1.1)\.
- \[17\]M\. Frank P\. Wolfeet al\.\(1956\)An algorithm for quadratic programming\.Naval research logistics quarterly3\(1\-2\),pp\. 95–110\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p3.1)\.
- \[18\]D\. Garber and O\. Meshi\(2016\)Linear\-memory and decomposition\-invariant linearly convergent conditional gradient algorithm for structured polytopes\.Advances in neural information processing systems29\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p3.1)\.
- \[19\]M\. X\. Goemans and D\. P\. Williamson\(1995\)Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming\.Journal of the ACM \(JACM\)42\(6\),pp\. 1115–1145\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p1.1)\.
- \[20\]K\. Helsgaun\(2017\)An extension of the lin\-kernighan\-helsgaun tsp solver for constrained traveling salesman and vehicle routing problems\.Roskilde: Roskilde University12,pp\. 966–980\.Cited by:[§D\.3](https://arxiv.org/html/2609.25728#A4.SS3.p3.1),[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p5.1)\.
- \[21\]M\. Jaggi\(2013\)Revisiting Frank–Wolfe: projection\-free sparse convex optimization\.InProceedings of the 30th International Conference on Machine Learning,Vol\.28,pp\. 427–435\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p3.1),[§4\.1](https://arxiv.org/html/2609.25728#S4.SS1.p4.1)\.
- \[22\]N\. Karalias and A\. Loukas\(2020\)Erdos goes neural: an unsupervised learning framework for combinatorial optimization on graphs\.Advances in Neural Information Processing Systems33,pp\. 6659–6672\.Cited by:[§B\.2](https://arxiv.org/html/2609.25728#A2.SS2.p4.1),[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[23\]N\. Karalias, A\. Rafiey, Y\. Xu, Z\. Luo, B\. Tahmasebi, C\. Jiang, and S\. Jegelka\(2025\)Geometric algorithms for neural combinatorial optimization with constraints\.InThe Thirty\-ninth Annual Conference on Neural Information Processing Systems,Cited by:[§B\.1](https://arxiv.org/html/2609.25728#A2.SS1.SSS0.Px1.p1.1),[§B\.1](https://arxiv.org/html/2609.25728#A2.SS1.SSS0.Px2.p1.1),[§B\.2](https://arxiv.org/html/2609.25728#A2.SS2.p4.1),[§B\.3](https://arxiv.org/html/2609.25728#A2.SS3.p1.1),[Appendix E](https://arxiv.org/html/2609.25728#A5.SS0.SSS0.Px1.p2.1),[Appendix G](https://arxiv.org/html/2609.25728#A7.SS0.SSS0.Px1.p2.1),[§1](https://arxiv.org/html/2609.25728#S1.p3.1),[§1](https://arxiv.org/html/2609.25728#S1.p5.1),[§1](https://arxiv.org/html/2609.25728#S1.p6.1),[§1](https://arxiv.org/html/2609.25728#S1.p7.1),[§2](https://arxiv.org/html/2609.25728#S2.p2.1),[§4\.2](https://arxiv.org/html/2609.25728#S4.SS2.p2.1),[§5\.1](https://arxiv.org/html/2609.25728#S5.SS1.p2.1)\.
- \[24\]N\. Karalias, J\. Robinson, A\. Loukas, and S\. Jegelka\(2022\)Neural set function extensions: learning with discrete functions in high dimensions\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Cited by:[§1](https://arxiv.org/html/2609.25728#S1.p3.1),[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[25\]A\. R\. Karlin, N\. Klein, and S\. O\. Gharan\(2021\)A \(slightly\) improved approximation algorithm for metric tsp\.InProceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing,pp\. 32–45\.Cited by:[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p2.1)\.
- \[26\]A\. R\. Karlin, N\. Klein, and S\. Oveis Gharan\(2023\)A deterministic better\-than\-3/2 approximation algorithm for metric tsp\.InInternational Conference on Integer Programming and Combinatorial Optimization \(IPCO\),pp\. 261–274\.Cited by:[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p2.1)\.
- \[27\]N\. Keriven and S\. Vaiter\(2023\)What functions can graph neural networks compute on random graphs? the role of positional encoding\.Advances in Neural Information Processing Systems \(NeurIPS\)36,pp\. 11823–11849\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[28\]E\. Khalil, H\. Dai, Y\. Zhang, B\. Dilkina, and L\. Song\(2017\)Learning combinatorial optimization algorithms over graphs\.InAdvances in Neural Information Processing Systems \(NeurIPS\),pp\. 6348–6358\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[29\]M\. Kim, J\. Park, and J\. Park\(2022\)Sym\-nco: leveraging symmetricity for neural combinatorial optimization\.Advances in Neural Information Processing Systems35,pp\. 1936–1949\.Cited by:[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p5.1)\.
- \[30\]V\. Konda and J\. Tsitsiklis\(1999\)Actor\-critic algorithms\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.12,pp\.\.Cited by:[§B\.2](https://arxiv.org/html/2609.25728#A2.SS2.p4.1)\.
- \[31\]W\. Kool, H\. van Hoof, and M\. Welling\(2019\)Attention, learn to solve routing problems\!\.InInternational Conference on Learning Representations \(ICLR\),Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[32\]R\. G\. Krishnan, S\. Lacoste\-Julien, and D\. Sontag\(2015\)Barrier frank\-wolfe for marginal inference\.Advances in Neural Information Processing Systems28\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p3.1)\.
- \[33\]Y\. Kushinsky, H\. Maron, N\. Dym, and Y\. Lipman\(2019\)Sinkhorn algorithm for lifted assignment problems\.SIAM Journal on Imaging Sciences12\(2\),pp\. 716–735\.Cited by:[§C\.3](https://arxiv.org/html/2609.25728#A3.SS3.p2.1)\.
- \[34\]S\. Lacoste\-Julien, M\. Jaggi, M\. Schmidt, and P\. Pletscher\(2013\)Block\-coordinate Frank\-Wolfe optimization for structural SVMs\.InProceedings of the 30th International Conference on Machine Learning \(ICML\),Proceedings of Machine Learning Research\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p3.1),[§4\.1](https://arxiv.org/html/2609.25728#S4.SS1.p4.1)\.
- \[35\]S\. Lacoste\-Julien and M\. Jaggi\(2015\)On the global linear convergence of frank\-wolfe optimization variants\.Advances in Neural Information Processing Systems \(NeurIPS\)28\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p3.1)\.
- \[36\]X\. Li, E\. Liang, and M\. Chen\(2025\)Gauge flow matching for efficient constrained generative modeling over general convex set\.InICLR 2025 Workshop on Deep Generative Model in Machine Learning: Theory, Principle and Efficacy,Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[37\]E\. Liang, M\. Chen, and S\. H\. Low\(2023\)Low complexity homeomorphic projection to ensure neural\-network solution feasibility for optimization over \(non\-\) convex set\.InInternational Conference on Machine Learning \(ICML\),pp\. 20623–20649\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[38\]E\. Liang, M\. Chen, and S\. H\. Low\(2024\)Homeomorphic projection to ensure neural\-network solution feasibility for constrained optimization\.Journal of Machine Learning Research25\(329\),pp\. 1–55\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[39\]E\. Liang and M\. ChenEfficient bisection projection to ensure nn solution feasibility for optimization over general set\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[40\]D\. Lim, J\. D\. Robinson, L\. Zhao, T\. Smidt, S\. Sra, H\. Maron, and S\. Jegelka\(2023\)Sign and basis invariant networks for spectral graph representation learning\.InThe Eleventh International Conference on Learning Representations \(ICLR\),Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[41\]C\. Liu, Z\. Jiang, R\. Wang, L\. Huang, P\. Lu, and J\. Yan\(2023\)Revocable deep reinforcement learning with affinity regularization for outlier\-robust graph matching\.InThe Eleventh International Conference on Learning Representations \(ICLR\),Cited by:[§C\.1](https://arxiv.org/html/2609.25728#A3.SS1.SSS0.Px6.p1.1),[§C\.2](https://arxiv.org/html/2609.25728#A3.SS2.p1.1),[§5\.2](https://arxiv.org/html/2609.25728#S5.SS2.p4.1)\.
- \[42\]F\. Locatello, R\. Khanna, M\. Tschannen, and M\. Jaggi\(2017\)A unified optimization view on generalized matching pursuit and frank\-wolfe\.InArtificial intelligence and statistics \(AISTAT\),pp\. 860–868\.Cited by:[Appendix A](https://arxiv.org/html/2609.25728#A1.p4.2.1),[§2](https://arxiv.org/html/2609.25728#S2.p3.1),[§4\.1](https://arxiv.org/html/2609.25728#S4.SS1.p4.1)\.
- \[43\]L\. Lovász\(1983\)Submodular functions and convexity\.InMathematical Programming The State of the Art: Bonn 1982,pp\. 235–257\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p1.1)\.
- \[44\]J\. Ma, W\. Pan, Y\. Li, and J\. Yan\(2025\)COExpander: adaptive solution expansion for combinatorial optimization\.InInternational Conference on Machine Learning,Cited by:[§D\.1](https://arxiv.org/html/2609.25728#A4.SS1.SSS0.Px1.p1.1),[§D\.2](https://arxiv.org/html/2609.25728#A4.SS2.SSS0.Px1.p1.1),[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p4.1)\.
- \[45\]Y\. Min, F\. Wenkel, M\. Perlmutter, and G\. Wolf\(2022\)Can hybrid geometric scattering networks help solve the maximum clique problem?\.Advances in Neural Information Processing Systems \(NeurIPS\)35,pp\. 22713–22724\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[46\]V\. Mirrokni, R\. P\. Leme, A\. Vladu, and S\. C\. Wong\(2017\)Tight bounds for approximate Carathéodory and beyond\.InProceedings of the 34th International Conference on Machine Learning \(ICML\),Vol\.70,pp\. 2440–2448\.Cited by:[§1](https://arxiv.org/html/2609.25728#S1.p4.1),[§2](https://arxiv.org/html/2609.25728#S2.p3.1),[§4\.1](https://arxiv.org/html/2609.25728#S4.SS1.p8.1)\.
- \[47\]K\. Mundinger, M\. Zimmer, A\. Kiem, C\. Spiegel, and S\. Pokutta\(2025\)Neural discovery in mathematics: do machines dream of colored planes?\.InForty\-second International Conference on Machine Learning \(ICML\),Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[48\]K\. Murota\(1998\)Discrete convex analysis\.Mathematical Programming83\(1\),pp\. 313–371\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p1.1)\.
- \[49\]G\. L\. Nemhauser, L\. A\. Wolsey, and M\. L\. Fisher\(1978\)An analysis of approximations for maximizing submodular set functions \- I\.Math\. Program\.14\(1\),pp\. 265–294\.Cited by:[§B\.2](https://arxiv.org/html/2609.25728#A2.SS2.p3.1)\.
- \[50\]R\. R\. Nerem, Z\. Luo, A\. Rafiey, and Y\. Wang\(2025\)Differentiable extensions with rounding guarantees for combinatorial optimization over permutations\.InThe Thirty\-ninth Annual Conference on Neural Information Processing Systems,Cited by:[Appendix G](https://arxiv.org/html/2609.25728#A7.SS0.SSS0.Px2.p1.2),[§1](https://arxiv.org/html/2609.25728#S1.p3.1),[§1](https://arxiv.org/html/2609.25728#S1.p7.1),[§2](https://arxiv.org/html/2609.25728#S2.p2.1),[§4\.2](https://arxiv.org/html/2609.25728#S4.SS2.p3.1)\.
- \[51\]G\. Pisier\(1981\)Remarques sur un résultat non publié de b\. maurey\.Séminaire d’Analyse fonctionnelle \(dit" Maurey\-Schwartz"\),pp\. 1–12\.Cited by:[§4\.1](https://arxiv.org/html/2609.25728#S4.SS1.p8.1)\.
- \[52\]R\. Qiu, Z\. Sun, and Y\. Yang\(2022\)Dimes: a differentiable meta solver for combinatorial optimization problems\.Advances in Neural Information Processing Systems \(NeurIPS\)35,pp\. 25531–25546\.Cited by:[§D\.3](https://arxiv.org/html/2609.25728#A4.SS3.p5.1),[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p5.1)\.
- \[53\]G\. Reinelt\(1991\)TSPLIB—a traveling salesman problem library\.ORSA journal on computing3\(4\),pp\. 376–384\.Cited by:[§D\.6](https://arxiv.org/html/2609.25728#A4.SS6.p2.1)\.
- \[54\]R\. Sato, M\. Yamada, and H\. Kashima\(2019\)Approximation ratios of graph neural networks for combinatorial problems\.Advances in Neural Information Processing Systems \(NeurIPS\)32\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[55\]M\. J\. Schuetz, J\. K\. Brubaker, and H\. G\. Katzgraber\(2022\)Combinatorial optimization with physics\-inspired graph neural networks\.Nature Machine Intelligence4\(4\),pp\. 367–377\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[56\]R\. Sinkhorn\(1964\)A relationship between arbitrary positive matrices and doubly stochastic matrices\.The annals of mathematical statistics35\(2\),pp\. 876–879\.Cited by:[§C\.1](https://arxiv.org/html/2609.25728#A3.SS1.SSS0.Px4.p2.1),[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[57\]H\. Sun, E\. K\. Guha, and H\. Dai\(2022\)Annealed training for combinatorial optimization on graphs\.arXiv preprint arXiv:2207\.11542\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[58\]Z\. Sun and Y\. Yang\(2023\)Difusco: graph\-based diffusion solvers for combinatorial optimization\.Advances in Neural Information Processing Systems \(NeurIPS\)36,pp\. 3706–3731\.Cited by:[§D\.3](https://arxiv.org/html/2609.25728#A4.SS3.p5.1),[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p5.1)\.
- \[59\]Z\. Tan and Y\. Mu\(2024\)Learning solution\-aware transformers for efficiently solving quadratic assignment problem\.InProceedings of the 41st International Conference on Machine Learning \(ICML\),pp\. 47627–47648\.Cited by:[§C\.1](https://arxiv.org/html/2609.25728#A3.SS1.SSS0.Px1.p1.1),[§C\.1](https://arxiv.org/html/2609.25728#A3.SS1.SSS0.Px6.p1.1),[§C\.2](https://arxiv.org/html/2609.25728#A3.SS2.p1.1),[§5\.2](https://arxiv.org/html/2609.25728#S5.SS2.p3.1),[§5\.2](https://arxiv.org/html/2609.25728#S5.SS2.p4.1)\.
- \[60\]M\. Tawarmalani and N\. V\. Sahinidis\(2002\)Convex extensions and envelopes of lower semi\-continuous functions\.Mathematical programming93\(2\),pp\. 247–263\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p1.1)\.
- \[61\]J\. Toenshoff, M\. Ritzert, H\. Wolf, and M\. Grohe\(2021\)Graph neural networks for maximum constraint satisfaction\.Frontiers in artificial intelligence3,pp\. 580607\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[62\]O\. Vinyals, M\. Fortunato, and N\. Jaitly\(2015\)Pointer networks\.InAdvances in Neural Information Processing Systems \(NeurIPS\),pp\. 2692–2700\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[63\]H\. P\. Wang, N\. Wu, H\. Yang, C\. Hao, and P\. Li\(2022\)Unsupervised learning for combinatorial optimization with principled objective relaxation\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[64\]R\. Wang, L\. Shen, Y\. Chen, X\. Yang, D\. Tao, and J\. Yan\(2022\)Towards one\-shot neural combinatorial solvers: theoretical and empirical notes on the cardinality\-constrained case\.InThe Eleventh International Conference on Learning Representations \(ICLR\),Cited by:[§B\.1](https://arxiv.org/html/2609.25728#A2.SS1.SSS0.Px1.p1.1),[§B\.2](https://arxiv.org/html/2609.25728#A2.SS2.p4.1),[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[65\]R\. Wang, J\. Yan, and X\. Yang\(2021\)Neural graph matching network: learning lawler’s quadratic assignment problem with extension to hypergraph and multiple\-graph matching\.IEEE Transactions on Pattern Analysis and Machine Intelligence44\(9\),pp\. 5261–5279\.Cited by:[§C\.2](https://arxiv.org/html/2609.25728#A3.SS2.p1.1)\.
- \[66\]R\. Wang, Y\. Zhang, Z\. Guo, T\. Chen, X\. Yang, and J\. Yan\(2023\)LinSATNet: the positive linear satisfiability neural networks\.InInternational Conference on Machine Learning \(ICML\),pp\. 36605–36625\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[67\]D\. P\. Williamson and D\. B\. Shmoys\(2011\)The design of approximation algorithms\.Cambridge university press\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p1.1)\.
- \[68\]Y\. Xia, X\. Yang, Z\. Liu, Z\. Liu, L\. Song, and J\. Bian\(2024\)Position: rethinking post\-hoc search\-based neural approaches for solving large\-scale traveling salesman problems\.InInternational Conference on Machine Learning \(ICML\),pp\. 54178–54190\.Cited by:[§D\.2](https://arxiv.org/html/2609.25728#A4.SS2.SSS0.Px3.p1.1),[Appendix E](https://arxiv.org/html/2609.25728#A5.SS0.SSS0.Px1.p2.1),[§5\.3](https://arxiv.org/html/2609.25728#S5.SS3.p4.1)\.
- \[69\]H\. Zeng, C\. Yang, Y\. Zhou, C\. Yang, and Q\. Guo\(2024\)GLinSAT: the general linear satisfiability neural network layer by accelerated gradient descent\.Advances in Neural Information Processing Systems \(NeurIPS\)37,pp\. 122584–122615\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p2.1)\.
- \[70\]B\. Zhou and Y\. Sun\(2022\)Approximate Frank–Wolfe algorithms over graph\-structured support sets\.InProceedings of the 39th International Conference on Machine Learning,Vol\.162,pp\. 27303–27337\.Cited by:[§2](https://arxiv.org/html/2609.25728#S2.p3.1),[§4\.1](https://arxiv.org/html/2609.25728#S4.SS1.p4.1)\.

## Appendix AProofs for[Theorem4\.3](https://arxiv.org/html/2609.25728#S4.Thmtheorem3)and Corollary[4\.4](https://arxiv.org/html/2609.25728#S4.Thmtheorem4)

###### Proof of[Theorem4\.3](https://arxiv.org/html/2609.25728#S4.Thmtheorem3)\.

The convex\-combination structure follows directly from the update𝐮t=\(1−γt\)​𝐮t−1\+γt​𝐯t\\mathbf\{u\}\_\{t\}=\(1\-\\gamma\_\{t\}\)\\mathbf\{u\}\_\{t\-1\}\+\\gamma\_\{t\}\\mathbf\{v\}\_\{t\}withγt∈\[0,1\]\\gamma\_\{t\}\\in\[0,1\]and the initialization𝐮0=𝐯0∈ext⁡\(𝒫\)\\mathbf\{u\}\_\{0\}=\\mathbf\{v\}\_\{0\}\\in\\operatorname\{ext\}\(\\mathcal\{P\}\)\. Inductively,

𝐮T−1=∑t=0T−1αt​𝐯t,𝐯t∈ext⁡\(𝒫\),αt≥0,∑t=0T−1αt=1\.\\mathbf\{u\}\_\{T\-1\}=\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}\\mathbf\{v\}\_\{t\},\\qquad\\mathbf\{v\}\_\{t\}\\in\\operatorname\{ext\}\(\\mathcal\{P\}\),\\qquad\\alpha\_\{t\}\\geq 0,\\qquad\\sum\_\{t=0\}^\{T\-1\}\\alpha\_\{t\}=1\.
For the convergence rate, letD=diam⁡\(𝒫\)D=\\mathrm\{diam\}\(\\mathcal\{P\}\)andJ⁡\(𝐮\)=12​‖𝐱−𝐮‖22J\(\\mathbf\{u\}\)=\\frac\{1\}\{2\}\\\|\\mathbf\{x\}\-\\mathbf\{u\}\\\|\_\{2\}^\{2\}\. The gradient ofJJis11\-Lipschitz, so the Frank–Wolfe curvature ofJJover𝒫\\mathcal\{P\}is bounded byD2D^\{2\}\. The standard Frank–Wolfe rate with exact linear minimization/maximization oracles and line search gives, afterkkupdates,

J⁡\(𝐮k\)−J⁡\(𝐮⋆\)≤2​D2k\+2,𝐮⋆∈arg​min𝐮∈𝒫⁡J​\(𝐮\)\.J\(\\mathbf\{u\}\_\{k\}\)\-J\(\\mathbf\{u\}^\{\\star\}\)\\leq\\frac\{2D^\{2\}\}\{k\+2\},\\qquad\\mathbf\{u\}^\{\\star\}\\in\\operatorname\*\{arg\\,min\}\_\{\\mathbf\{u\}\\in\\mathcal\{P\}\}J\(\\mathbf\{u\}\)\.Settingk=T−1k=T\-1yields

J⁡\(𝐮T−1\)−J⁡\(𝐮⋆\)≤2​D2T\+1\.J\(\\mathbf\{u\}\_\{T\-1\}\)\-J\(\\mathbf\{u\}^\{\\star\}\)\\leq\\frac\{2D^\{2\}\}\{T\+1\}\.
It remains to justify the almost\-everywhere differentiability claim under the stable oracle decision assumption\. For almost every𝐱\\mathbf\{x\}, the assumption gives a neighborhood on which all selected vertices and line\-search clipping regimes are fixed\. Fix such an input and neighborhood\. The initialization𝐮0=𝐯0\\mathbf\{u\}\_\{0\}=\\mathbf\{v\}\_\{0\}is constant on this neighborhood\. Recursively, if𝐮t−1​\(𝐱\)\\mathbf\{u\}\_\{t\-1\}\(\\mathbf\{x\}\)is smooth there, then

𝐫t​\(𝐱\)=𝐱−𝐮t−1​\(𝐱\),𝐝t​\(𝐱\)=𝐯t−𝐮t−1​\(𝐱\)\\mathbf\{r\}\_\{t\}\(\\mathbf\{x\}\)=\\mathbf\{x\}\-\\mathbf\{u\}\_\{t\-1\}\(\\mathbf\{x\}\),\\qquad\\mathbf\{d\}\_\{t\}\(\\mathbf\{x\}\)=\\mathbf\{v\}\_\{t\}\-\\mathbf\{u\}\_\{t\-1\}\(\\mathbf\{x\}\)are smooth because𝐯t\\mathbf\{v\}\_\{t\}is fixed\. When𝐝t​\(𝐱\)≠0\\mathbf\{d\}\_\{t\}\(\\mathbf\{x\}\)\\neq 0, the unclipped line\-search value

γ~t​\(𝐱\)=⟨𝐫t​\(𝐱\),𝐝t​\(𝐱\)⟩‖𝐝t​\(𝐱\)‖22\\tilde\{\\gamma\}\_\{t\}\(\\mathbf\{x\}\)=\\frac\{\\langle\\mathbf\{r\}\_\{t\}\(\\mathbf\{x\}\),\\mathbf\{d\}\_\{t\}\(\\mathbf\{x\}\)\\rangle\}\{\\\|\\mathbf\{d\}\_\{t\}\(\\mathbf\{x\}\)\\\|\_\{2\}^\{2\}\}is a smooth algebraic function on a possibly smaller neighborhood with nonzero denominator\. When𝐝t=0\\mathbf\{d\}\_\{t\}=0, we use the algorithmic conventionγt=0\\gamma\_\{t\}=0, and the update is locally constant in that direction\. Since the clipping regime is fixed on the neighborhood,γt\\gamma\_\{t\}is either00,11, orγ~t\\tilde\{\\gamma\}\_\{t\}, and is therefore differentiable there\. Hence

𝐮t=\(1−γt\)​𝐮t−1\+γt​𝐯t\\mathbf\{u\}\_\{t\}=\(1\-\\gamma\_\{t\}\)\\mathbf\{u\}\_\{t\-1\}\+\\gamma\_\{t\}\\mathbf\{v\}\_\{t\}involves only differentiable algebraic operations\. The same recursive update shows that each coefficientαt​\(𝐱\)\\alpha\_\{t\}\(\\mathbf\{x\}\)is differentiable on the neighborhood\. The only excluded inputs lie in the zero\-volume switching set from the stable oracle decision assumption, so the coefficient map𝐱↦𝜶⁡\(𝐱\)\\mathbf\{x\}\\mapsto\\bm\{\\alpha\}\(\\mathbf\{x\}\)is differentiable almost everywhere\. ∎

###### Proof of Corollary[4\.4](https://arxiv.org/html/2609.25728#S4.Thmtheorem4)\.

The convex\-combination conclusion is unchanged because the update still usesγt∈\[0,1\]\\gamma\_\{t\}\\in\[0,1\]and vertices of𝒫\\mathcal\{P\}\. The relative FW\-gap condition

⟨𝐫t,𝐯t−𝐮t−1⟩≥δ​max𝐯∈𝒫​⟨𝐫t,𝐯−𝐮t−1⟩\\langle\\mathbf\{r\}\_\{t\},\\mathbf\{v\}\_\{t\}\-\\mathbf\{u\}\_\{t\-1\}\\rangle\\geq\\delta\\max\_\{\\mathbf\{v\}\\in\\mathcal\{P\}\}\\langle\\mathbf\{r\}\_\{t\},\\mathbf\{v\}\-\\mathbf\{u\}\_\{t\-1\}\\rangleis the maximization form of the standard approximate Frank–Wolfe oracle condition\. Applying the approximate Frank–Wolfe theorem of[Locatello et al\. \[42\]](https://arxiv.org/html/2609.25728#bib.bib51)toJ⁡\(𝐮\)=12​‖𝐱−𝐮‖22J\(\\mathbf\{u\}\)=\\frac\{1\}\{2\}\\\|\\mathbf\{x\}\-\\mathbf\{u\}\\\|\_\{2\}^\{2\}, with Lipschitz constant11and diameterDD, gives afterT−1T\-1updates

J⁡\(𝐮T−1\)−J⁡\(𝐮⋆\)≤2​\(D2/δ\+h0\)δ⁡\(T−1\)\+2,h0=J⁡\(𝐮0\)−J⁡\(𝐮⋆\)\.J\(\\mathbf\{u\}\_\{T\-1\}\)\-J\(\\mathbf\{u\}^\{\\star\}\)\\leq\\frac\{2\\left\(D^\{2\}/\\delta\+h\_\{0\}\\right\)\}\{\\delta\(T\-1\)\+2\},\\qquad h\_\{0\}=J\(\\mathbf\{u\}\_\{0\}\)\-J\(\\mathbf\{u\}^\{\\star\}\)\.The almost\-everywhere differentiability statement follows from the same stable oracle decision argument used above\. ∎

## Appendix BDetails and Further Experiments for Maximum Coverage

#### Hardware\.

All experiments, including training and inference on Maximum Coverage, TSP and QAP, are done on 16 cores \(32 threads\) of Intel\(R\) Xeon\(R\) Platinum 8268 CPU \(24 cores, 48 threads in total\), 32 GB DDR4 ram, with a single Nvidia RTX8000 48GB GPU\.

### B\.1Experiment Setup

#### Datasets\.

Following prior work\[[7](https://arxiv.org/html/2609.25728#bib.bib39),[64](https://arxiv.org/html/2609.25728#bib.bib50),[23](https://arxiv.org/html/2609.25728#bib.bib11)\], we evaluate on synthetic and real\-world bipartite graphs\(U,V,E\)\(U,V,E\), whereVVis the ground set and the task is to selectkknodes fromVVto maximize the total weight of covered nodes inUU\. For synthetic data, we use the Random Uniform dataset under two scales,Random500andRandom1000\. Each setting contains a training dataset of100100independently generated instances with\|V\|=500,\|U\|=1000\|V\|=500,\|U\|=1000\(Random500\) or\|V\|=1000,\|U\|=2000\|V\|=1000,\|U\|=2000\(Random1000\)\. Eachu∈Uu\\in Uis assigned an integer weight sampled uniformly from\{1,…,100\}\\\{1,\\ldots,100\\\}, and eachv∈Vv\\in Vcovers a random subset ofUUwhose cardinality is sampled uniformly from\{10,…,30\}\\\{10,\\ldots,30\\\}\. Testing are conducted on datasets sampled from the same distribution\. For real\-world testing, we additionally evaluate on three*Railway*instances derived from Italian railway crew assignment graphs:rail507\(\|V\|=507,\|U\|=63009\|V\|=507,\|U\|=63009\),rail516\(\|V\|=516,\|U\|=47311\|V\|=516,\|U\|=47311\), andrail582\(\|V\|=582,\|U\|=55515\|V\|=582,\|U\|=55515\)\.

#### Training\.

In general, we use the same training setup as\[[23](https://arxiv.org/html/2609.25728#bib.bib11)\], using a 3\-layer GNN with residual connection\. To enhance optimization stability and solution quality, we incorporate specific regularization strategies to promote exploration within the hypersimplex\. An entropy regularization term is applied with a cosine\-decayed weightλt\\lambda\_\{t\}, which is initialized at0\.050\.05and annealed to zero by the 30th epoch\. A dropout rate of0\.10\.1is also applied exclusively during the training phase\.

The model is trained for a total of8080epochs using the AdamW optimizer with a weight decay of1×10−41\\times 10^\{\-4\}\. We employ a batch size of44and a fixed random seed of4242to ensure reproducibility\. The learning rate follows awarmup\_cosinescheduler, which includes5050warmup epochs and decays from a peak of5×10−35\\times 10^\{\-3\}to a minimum of5×10−55\\times 10^\{\-5\}\. Data loading and preprocessing are handled efficiently by1616CPU worker threads\.

#### Inference\.

In inference, we run ourFW Decompositionfor 50 iterations for each instance, and select the best result out of the decomposed sets\. We do not use any randomized or local improvment method\.

### B\.2Baselines

We compare against three classes of baselines that span exact solvers, classical heuristics, and neural approaches\.

\(i\) Exact solvers \(non\-learning\)\.We employGurobito obtain exact MIP\-based solutions, imposing a time limit of 120 seconds for each graph\.

\(ii\) Classical heuristics \(non\-learning\)\.We include aRandombaseline, which samples subsets of sizekkuniformly at random over multiple trials and selects the best solution found within 240 seconds\. We also report theGreedy algorithm\[[49](https://arxiv.org/html/2609.25728#bib.bib60)\], which iteratively adds the element with the largest marginal gain to achieve the well\-known\(1−1e\)\(1\-\\tfrac\{1\}\{\\mathrm\{e\}\}\)\-approximation guarantee\. Additionally, we classifyUCOM2\[[7](https://arxiv.org/html/2609.25728#bib.bib39)\]as a non\-learning baseline due to its use of incremental greedy de\-randomization, and we report results for its short variant\.

\(iii\) Neural baselines\.We evaluate a diverse set of neural methods\. We includeCardNN\[[64](https://arxiv.org/html/2609.25728#bib.bib50)\]variants \(without Test\-Time Optimization, following\[[7](https://arxiv.org/html/2609.25728#bib.bib39)\]\) and a naive version ofEGN\[[22](https://arxiv.org/html/2609.25728#bib.bib14)\]with naive probabilistic objective construction and iterative rounding\. We also compare againstGeoNCO\[[23](https://arxiv.org/html/2609.25728#bib.bib11)\], which employs the same Carathéodory decomposition framework as our method but lacks our proposed enhancements\. Finally, we include a Reinforcement Learning \(RL\) baseline following\[[23](https://arxiv.org/html/2609.25728#bib.bib11)\]trained with the Actor\-Critic\[[30](https://arxiv.org/html/2609.25728#bib.bib59)\]algorithm on the same problem instances\.

### B\.3Ablation study on projection

We ablate the effect of explicitly projecting the network output back into the feasible solution polytope before runningFW Decomposition\.We use the same projection method as\[[23](https://arxiv.org/html/2609.25728#bib.bib11)\], while they rely on such a projection to enforce feasibility, our framework does not require it since feasibility is already guaranteed by the decomposition\.

Table[4](https://arxiv.org/html/2609.25728#A2.T4)shows that removing projection preserves solution quality and in fact yields slightly better coverage across bothrand500andrand1000settings, with slightly lower wall\-clock inference time\. We attribute this to the fact that the projection is a hard, piecewise operator that can compress the effective search region and reduce gradient signal diversity; in our pipeline, it becomes largely redundant because feasibility is enforced by the decomposition itself\. By contrast, allowing unconstrained continuous scores retains a richer exploration space for the model, whileFW Decompositionstill produces feasible candidates, leading to steadier optimization and marginally improved performance in practice\.

Table 4:Ablation study results\. Grouped by dataset size\.
### B\.4Ablation on Hard Instances for Greedy

#### Hard example for greedy\.

The greedy algorithm performs well, especially on random graphs, for the MC problem because the coverage function is monotone and submodular\. To show our method works even when greedy fails, we construct theGreedyTrap500andGreedyTrap1000datasets for MC\. In this construction, the greedy algorithm is pushed near its theoretical1−1/e=0\.632121\-1/e=0\.63212lower bound and attracted by sets with large initial marginal coverage, falling significantly short of the optimal solution\.

The construction is as follows\. We first choose a coverable subset of elements from𝒰\\mathcal\{U\}and partition it intokkdisjoint parts\. For each part, we add one set that covers the entire part; selecting thesekksets covers all coverable elements\. We then addkkgreedy\-trap sets\. Each trap set covers slightly more than a1/k1/kfraction of every part, so these sets have large marginal gain at the beginning and are selected by greedy\. The remaining sets are low\-degree filler sets\. ForGreedyTrap500, we use 500 sets and 1000 elements; forGreedyTrap1000, we use 1000 sets and 2000 elements\. In testing, we report the coverage ratio over all elements\.

Table 5:Performance comparison on theGreedyTrapdatasets\. Values represent coverage ratios in percentage averaged across 64 test instances\.As shown in Table[5](https://arxiv.org/html/2609.25728#A2.T5), greedy falls significantly short of OPT on these hard instances, whereas our method consistently outperforms greedy and achieves coverage much closer to OPT\. This shows that our method can still learn useful problem structure and maintain strong performance on data distributions where the standard greedy heuristic fails\.

### B\.5Effect of the Frank–Wolfe Iteration Budget

We investigate the effect of the number of Frank–Wolfe decomposition iterationsTTon solution quality and runtime for Maximum Coverage\. Each iteration of Algorithm[1](https://arxiv.org/html/2609.25728#alg1)performs one LMO call and a closed\-form line\-search step\. Note that the runtime is linear inTT\.

Table[6](https://arxiv.org/html/2609.25728#A2.T6)reports results for four configurations, with solution quality normalized by the weighted coverage obtained atT=200T=200\. IncreasingTTconsistently improves quality, with diminishing gains at larger budgets\. AtT=50T=50, the method already achieves 99\.391–99\.652% of the reference quality across all four configurations\. Increasing the budget toT=200T=200improves relative quality by only 0\.348–0\.609 percentage points, while increasing the reported runtime by approximately3\.73\.7–3\.8×3\.8\\times\. These results support our choice ofT=50T=50as a practical balance between solution quality and computational cost\.

Table 6:Maximum coverage quality and inference time under different FW iteration budgets\. Quality is normalized by the weighted coverage obtained atT=200T=200and reported as a percentage\. Time recorded on an NVIDIA GeForce RTX 4090 GPU and not directly comparable with other results\.

## Appendix CDetails and Further Experiments on QAP

In this section, we provide a detailed discussion of the numerical results for QAP\. The full experimental results are reported in[Table7](https://arxiv.org/html/2609.25728#A3.T7)\.

Table 7:Complete experiment results on QAPLIB\. The mean/max/min gaps are reported for each category\. The mean performance over all categories and the inference time \(s\) per instance are reported\.### C\.1Experiment Setup

#### Instance distribution and data splits\.

For synthetic training, we follow the data generation method of\[[59](https://arxiv.org/html/2609.25728#bib.bib52)\]and generate a synthetic training pool with 5120 instances,n=32n=32andp=0\.3p=0\.3\. During training, each gradient step samples a mini\-batch uniformly from the training pool\.

For test, we use the QAPLIB dataset to demostrate the generalization ability of our method\. QAPLIB\[[8](https://arxiv.org/html/2609.25728#bib.bib61)\]is a widely used benchmark for QAP, collecting real\-world and synthetic instances from multiple sources\. Each instance is specified by a flow matrix and a distance matrix, and the naming convention uses a class prefix \(e\.g\.,bur,chr,tai\) together with the problem size\. The distributions and structural properties of the flow/distance matrices vary substantially across classes, making QAPLIB a challenging testbed for cross\-class generalization\.

#### Input representation\.

GivenAAandBB, we build a node embedding for each matrix separately using per\-node statistics \(diagonal entry and row/column aggregations\) together with random\-projection features computed with a fixed Gaussian projection matrixRR\. Features are computed after per\-instance magnitude normalization \(i\.e\., dividing bymaxi,j⁡\|Mi​j\|\\max\_\{i,j\}\|M\_\{ij\}\|\) to stabilize training across varying scales\.

#### Model, optimizer, and training schedule\.

We use a two\-layer GraphSAGE aggregator with a residual skip connection and hidden sizemodel\_dim=128, applied independently toAAandBB\. The neural network outputs a matrix followed by a Sinkhorn projection layer\. We train with AdamW \(learning rate3⋅10−43\\cdot 10^\{\-4\}, weight decay10−410^\{\-4\}\) and gradient clipping at norm1\.01\.0\. Training runs for 100 epochs; each epoch consists of 100 instances randomly sampled from the training pool\. A single global RNG seed \(0\) is used for both training and inference\.

#### Sinkhorn as a practical projection layer\.

Although our FW\-based pipeline does*not*theoretically require the network output to lie inside the Birkhoff polytope \(theFW Decompositioncan be applied to any real matrix\), in QAP this relaxation is empirically fragile because QAP instances are substantially more heterogeneous than the other tasks: QAPLIB in particular exhibits highly dispersed scales, sparsity patterns, and instance\-generation mechanisms across data classes, leading to large distribution gaps relative to synthetic pretraining\. In this regime, leaving the model output unconstrained often yields score matrices whose magnitude and row/column mass distribution are poorly calibrated, which in turn destabilizes the downstream greedy steps and makes the distribution returned byFW Decompositionless informative\.

To mitigate this, we insert a Sinkhorn\[[56](https://arxiv.org/html/2609.25728#bib.bib49)\]projection layer with 20 iterations as a soft projection that maps raw assignment logits to an approximately doubly\-stochastic matrix\. This operation enforces nonnegativity and approximately unit row/column sums, placing𝐗\\mathbf\{X\}, the output of neural network, closer to the Birkhoff polytope\. It serves as a lightweight output\-range calibration mechanism that significantly improves numerical stability under the large inter\-class distribution shift in QAPLIB\. We emphasize that this Sinkhorn step is used for practical robustness, not as a theoretical requirement of theFW Decompositionitself\.

#### Decomposition settings\.

We run ourFW Decompositionalgorithm over permutation matrices with a iteration budgetT=2000T=2000in training andT=3000T=3000in inference\.

At each iteration, the LMO is instantiated by a*greedy*maximum\-weight assignment routine \(no Hungarian\), applied to the current residual matrix\. The resulting decomposition yields a sequence of permutations and their weights; in inference, we select the minimum\-cost permutation under the QAP objective for reporting\.

Table 8:QAPLIB results for instances withn\>64n\>64\.Instancepred\_bestbest\_knowngaplipa70a1729731697551\.90%lipa70b5872840460320027\.58%sko72759526625614\.63%lipa80a2577692531951\.81%lipa80b9983972776396228\.59%tai80b105783622481841504329\.25%tai80a153116641349918413\.43%sko811042149099814\.52%lipa90a3663293606301\.58%lipa90b160171971249044128\.24%sko9013207811553414\.32%sko100a17404415200214\.50%sko100b17437815389013\.31%sko100c17027614786215\.16%sko100d17057414957614\.04%sko100e17074614915014\.48%sko100f16955814903613\.77%tai100a236380542105246612\.28%tai100b1498196608118599613726\.32%wil1002961362730388\.46%esc12824464281\.25%tai150b61955628849889664324\.19%tho1509579958813339817\.79%tai256c507530964475929413\.39%Average gap26\.87%
#### QAPLIB evaluation\.

Following previous approaches\[[59](https://arxiv.org/html/2609.25728#bib.bib52),[41](https://arxiv.org/html/2609.25728#bib.bib53)\], We report per\-class gap statistics and the aggregate summary as in[Table7](https://arxiv.org/html/2609.25728#A3.T7)for instances of size 12\-64 since many baseline methods will run out of time on very large instances\. However, our method is also capable of solving instances with larger sizes, as listed below in[Table8](https://arxiv.org/html/2609.25728#A3.T8)\.

### C\.2Related Work \(QAP / Matching\)\.

QAP and matching are among the most challenging CO problems for learning\-based methods, largely due to severe distribution shifts and outliers such as those in QAPLIB\.NGM\[[65](https://arxiv.org/html/2609.25728#bib.bib62)\]was one of the earliest attempts to tackle QAP/graph matching with neural models\. Subsequent work proposedRGM\[[41](https://arxiv.org/html/2609.25728#bib.bib53)\], which introduced a*learn\-to\-construct \(L2C\)*pipeline for QAP, building a permutation by making a sequence of constructive decisions, following that,SAWT\[[59](https://arxiv.org/html/2609.25728#bib.bib52)\]introduced a*learn\-to\-improve \(L2I\)*pipeline, which starts from an initial permutation and repeatedly applies local modifications to improve it\.

Both L2C and L2I are inherently iterative: they start from an empty/partial solution or an initial solution and improve it step by step\. This makes inference naturally less efficient and less friendly to large instances, and robustness to outliers is achieved through relatively complex neural architectures and nontrivial inference procedures\.

Our general framework applies to QAP in a single\-pass direct\-output manner: we produce assignment logits once and deterministically output a permutation without repeated optimization loops, yielding an immediate speed advantage\. Moreover, with only a lightweight Sinkhorn normalization layer added to the base framework, we obtain outlier\-friendly behavior at very low additional cost\.

### C\.3Baselines

We compare against two classes of baselines that cover both classical matching solvers and representative neural methods on QAPLIB\.

\(i\) Classical solvers \(non\-learning\)\.We includeSM,RRWM, andSK\-JA\[[33](https://arxiv.org/html/2609.25728#bib.bib63)\]as standard non\-learning baselines for graph matching / QAP\.

\(ii\) Neural methods\.We evaluateNGM,RGM, andSAWTas representative learning\-based baselines on QAPLIB\.

We follow the same per\-class reporting protocol as prior QAPLIB evaluations, and report the mean/min/max optimality gaps for each QAPLIB class together with the average per\-instance inference time\. As the public releases of bothSAWTandRGMdo not include a complete evaluation setup and code for QAPLIB, we are unable to reproduce their experiments locally\. Instead, we directly cite their reported timing results\. Consequently, due to hardware disparities, these runtime comparisons should be interpreted as indicative references rather than strictly controlled benchmarks\. We also note an additional caveat in their QAPLIB results: for the WIL class they report multiple results, whereas within range \(n∈\[12,64\]n\\in\[12,64\]\) where they reported, WIL contains only a single instance \(WIL50\)\. Since we cannot verify the exact evaluation protocol \(e\.g\., whether additional variants or sizes were included\), we therefore quote their reported numbers in full without further reinterpretation\.

### C\.4Effect of the Frank–Wolfe Iteration Budget

We evaluate the sensitivity of QAP solution quality and runtime to the Frank–Wolfe decomposition budget \(TT\) on QAPLIB\. Table[9](https://arxiv.org/html/2609.25728#A3.T9)reports relative quality normalized by the weighted solution obtained atT=10000T=10000, together with runtime for each budget\.

Solution quality improves substantially as the budget increases fromT=50T=50toT=1000T=1000, rising from 94\.247% to 99\.237% of the reference quality\. Further iterations yield progressively smaller improvements\. At our chosen budget ofT=3000T=3000, relative quality reaches 99\.837%, with a runtime of 4\.410 seconds\. Increasing the budget toT=10000T=10000adds only 0\.163 percentage points of relative quality while increasing runtime to 15\.798 seconds, approximately3\.6×3\.6\\timeshigher\. We therefore useT=3000T=3000to balance solution quality and runtime\.

Compared with Maximum Coverage, QAP empirically benefits from a substantially larger decomposition budget\. The approximation guarantee in Theorem[4\.3](https://arxiv.org/html/2609.25728#S4.Thmtheorem3)provides qualitative context: the squared approximation error is bounded by4​D2/\(T\+1\)4D^\{2\}/\(T\+1\), whereDDis the diameter of the feasible polytope\. Our choice ofTTis guided by this and the observed quality–runtime trade\-off\.

Table 9:QAP solution quality under different FW iteration budgets on QAPLIB\. Results are normalized by the weighted solution obtained at T = 10000\. Time recorded on an NVIDIA GeForce RTX 4090 GPU and not directly comparable with other results\.

## Appendix DDetails and Further Experiments on TSP

### D\.1Experiment Setup

#### Instance distribution and data splits\.

We train and evaluate on the 2D Euclidean TSP, where each instance consists ofnnpoints in\[0,1\]2\[0,1\]^\{2\}and tour length is computed under standard Euclidean distances\. We use the training and evaluation data fromCOExpander\[[44](https://arxiv.org/html/2609.25728#bib.bib19)\]\. During training, for a givennn, at each epoch, 64 instances are sampled from the training pool, and different data are used across epochs\. The test split has 1280 instances for TSP\-50 and TSP\-100, 128 instances for TSP\-500 and 32 instances for TSP\-1000, so that most methods could complete the evaluation in reasonable time frame\.

#### Graph construction and input features\.

For each instance we construct the complete graph onnnnodes \(undirected edges enumerated withu<vu<v\), and represent it in PyG using the standard bidirected encoding \(each undirected edge appears in both directions inedge\_index\)\. Node features are the raw coordinates augmented with a Fourier positional encoding \(withF=4F\{=\}4frequencies\) passed through a small MLP before message passing\. Edge features consist of a single scalar channel: the LP relaxation valuexLP​\(e\)x\_\{\\text\{LP\}\}\(e\)for the undirected edge, duplicated across both directed copies\.

#### LP relaxation feature computation\.

We computexLPx\_\{\\text\{LP\}\}by solving the subtour\-elimination LP using the QSopt LP solver, which is also the LP solver used by Concorde\. Edge costs are obtained from Euclidean distances and scaled to integer weights using a fixed coordinate scale of10610^\{6\}for the QSopt interface\. LP solves are performed on CPU for every instance in both train and test splits\.

#### Model, optimizer, and training schedule\.

The encoder is a GATv2 model with residual blocks: hidden size 128, 3 message\-passing layers, and 4 attention heads per layer\. We train with AdamW with a learning rate of10−310^\{\-3\}and a weight decay of10−410^\{\-4\}for 50 epochs, using a single global RNG seed \(0\) for reproducibility\. The best model is saved and used for inference\.

#### Decomposition hyperparameters used in training and inference\.

We use a fixed decomposition budget ofT=64T=64iterations in both training and testing, and a max\-entropy regularization with temperaturefw\_entropy\_tau=2\.0=2\.0is added to the decomposition\. All reported TSP test results use the same hyperparameters as training and the reported result is the best permutation from the decomposition\.

### D\.2Related Work: Neural Methods for TSP

#### Learning paradigms: LC vs\. GP\.

Neural approaches for TSP are often categorized by the granularity at which they make decisions\. According to\[[44](https://arxiv.org/html/2609.25728#bib.bib19)\], the*Learning\-to\-Construct \(LC\)*paradigm builds a tour sequentially, repeatedly selecting the next decision conditioned on the current partial solution\. A key strength of LC is that feasibility can be*hard\-enforced*during decoding \(e\.g\., by imposed masks that rule out invalid moves\), which provides a transparent and reliable path to valid tours\. However, this same sequential structure typically incurs high inference latency on large instances, since the model must execute many dependent steps\.

In contrast, the*Global Prediction \(GP\)*paradigm predicts a globally complete structure in one \(or a few\) forward passes, commonly as an edge\-probability heatmap, and then applies a separate recovery procedure to obtain an integer\-feasible tour\. GP can be substantially faster at inference and has become a dominant template for large\-scale settings, including diffusion\- and consistency\-based solvers that generate global solutions through iterative denoising\. At the same time, GP shifts difficulty into the heatmap\-to\-solution stage: the final performance can depend heavily on the specific decoding/recovery algorithm \(e\.g\., Parallel Sampling, Greedy, Monte Carlo Tree Search \(MCTS\)\) and additional improvement heuristics \(e\.g\., 2\-opt\), making the pipeline sensitive to post\-processing choices and sometimes blurring what is being evaluated \(the neural predictor vs\. the downstream solver\), since heatmap recovery is itself a combinatorial optimization subroutine\.

#### Adaptive Expansion \(AE\) and the remaining heatmap bottleneck\.

COExpanderintroduces*Adaptive Expansion \(AE\)*as an intermediate paradigm intended to combine the best of both worlds: it uses a global predictor to produce a heatmap under partial\-solution prompting, and then expands the determined set via a determination/decoding operator, effectively adapting the decision granularity across iterations\. While AE improves the efficiency–feasibility trade\-off relative to pure LC or pure GP, it still inherits a core limitation of GP\-style pipelines: it relies on predicting a heatmap and then recovering an integer tour via post\-processing\. So overall success remains coupled to the robustness and tuning of heatmap\-based recovery\.

#### Toward end\-to\-end pipelines with reduced problem\-specific setup\.

[Xia et al\. \[68\]](https://arxiv.org/html/2609.25728#bib.bib64)critically examine the prevalent*heatmap\-guided post\-hoc search*paradigm for large\-scale TSP \(notably heatmap\-guided MCTS\), questioning the practical value of learning heatmaps that ultimately require substantial downstream search and advocating a shift toward more*autonomous and generalizable*ML pipelines with less hand\-crafted machinery and stronger guarantees\. Motivated by this perspective, our method aims to retain LC’s feasibility discipline and GP’s efficiency while avoiding the heatmap bottleneck: we adopt an end\-to\-end pipeline that differentiates through the selected FW branch, with the*optimization objective aligned with the prediction target*—we do not learn a surrogate heatmap and then rely on a separate recovery/search step, but instead directly optimize for improved feasible tours through a structured, deterministic procedure\. This design reduces reliance on hand\-crafted post\-processing, improves robustness across instance distributions, and enables stronger theoretical characterization of the resulting algorithmic guarantees\.

### D\.3Baselines

We compare against four classes of baselines that span both classical optimization and representative neural paradigms under different supervision regimes\.

\(i\) Exact solvers \(non\-learning\)\.We reportConcordeas an exact TSP solver, and theGurobioptimizer as an exact baseline on small instances\. For larger instances, exact results are omitted when the solver fails to finish within the runtime cutoff\.

\(ii\) Classical heuristics \(non\-learning\)\.We includeLKH3\[[20](https://arxiv.org/html/2609.25728#bib.bib66)\]as a strong hand\-engineered heuristic, run for 500 trials\.

\(iii\) Neural LC baselines\.We evaluateRL4CO \(SymNCO\)\[[5](https://arxiv.org/html/2609.25728#bib.bib67)\]as a representative*Learning\-to\-Construct \(LC\)*method trained with reinforcement learning\. LC\-style solvers typically decode solutions sequentially and rely on imposed masks \(or equivalent constraints\) to hard\-guarantee feasibility during construction\.

\(iv\) Neural GP/AE baselines\.We test heatmap\-driven pipelines in both RL\-style and supervised regimes\.DIMES\[[52](https://arxiv.org/html/2609.25728#bib.bib58)\]is included as a representative GP\-style method trained via meta\-RL\. On the supervised side, we includeDIFUSCO\[[58](https://arxiv.org/html/2609.25728#bib.bib57)\], a diffusion\-based GP approach, andCOExpander, a heatmap\-based*Adaptive Expansion \(AE\)*method that achieves state\-of\-the\-art performance\.

Overall, under a fixed compute budget, we aim to test the landscape as comprehensively as possible across \(a\) paradigms \(LC vs\. GP vs\. AE\) and \(b\) supervision regimes \(RL / meta\-RL vs\. supervised learning\), while transparently reporting any omissions caused by solver timeouts or training\-cost constraints\.

### D\.4Detailed Results

Table 10:Raw performance comparison on TSP instances\. We report the objective value \(Obj\.\) and inference time \(s\)\.∗Results on larger instances are omitted since it was unable to finish within the cutoff time\. †Results are omitted due to prohibitive computational costs during training on large datasets\.

#### Runtime breakdown\.

In[Table11](https://arxiv.org/html/2609.25728#A4.T11), we measure the average per\-instance wall\-clock cost of each stage in the TSP pipeline \. The Subtour\-Elimination LP andFW Decompositiondominates the runtime; in the learned matching version, the matching takes up a large portion of the overall runtime\.

Table 11:Average per\-instance time on TSP1000\.
#### Generalization Experiments \(Across Instance Sizes\)\.

We further evaluate how well a model trained on one problem size transfers to other sizes\. Specifically, we train three separate models on TSP\-100, TSP\-500, and TSP\-1000, respectively, and test each model on the same evaluation splits at sizes\{100,500,1000\}\\\{100,500,1000\\\}\. Table[12](https://arxiv.org/html/2609.25728#A4.T12)reports the resulting tour lengths\.

Overall, we observe strong cross\-size generalization: performance is largely stable when transferring between sizes, with only modest degradation when training on smaller instances and testing on larger ones\. For example, a model trained on TSP\-100 remains competitive on TSP\-500/1000, albeit with a small gap relative to models trained on the target size\. Conversely, models trained on larger instances transfer well to smaller instances, with only minor differences on TSP\-100\. These results suggest that the learned scoring function captures size\-agnostic geometric regularities of Euclidean TSP, enabling effective transfer across instance scales without retraining on each target size\.

Table 12:Generalization performance across different instance sizes\.

### D\.5Discussion on Learned Matching

#### Motivation and positioning\.

In our TSP pipeline, each candidate spanning tree produced by the FW decomposition is rounded into a tour via the Christofides procedure\. A key step in Christofides is computing a perfect matching on the odd\-degree vertex set induced by the tree\. Our default implementation uses a greedy matching routine for efficiency, while the learned\-matching variant replaces this step with a trainable module that is*optimized end\-to\-end for the TSP objective*, rather than supervised to recover an exact minimum\-weight perfect matching\.

#### Matching as another structured polytope component\.

For a given tree, letOObe its odd\-degree vertex set \(with\|O\|\|O\|even\), and letEOE\_\{O\}denote all unordered pairs inOO, which forms a structured polytope\. The matching network outputs a vectorS∈\(0,1\)\|EO\|S\\in\(0,1\)^\{\|E\_\{O\}\|\}overEOE\_\{O\}\. We then run ourFW Decompositionprocedure on the matching polytope overOOto obtain a short convex decompositionAmatch≈∑kβk​MkA\_\{\\text\{match\}\}\\approx\\sum\_\{k\}\\beta\_\{k\}\\,\{M\_\{k\}\}, where eachMkM\_\{k\}is an*integral*perfect matching \(an extreme point of the polytope\),βk≥0\\beta\_\{k\}\\geq 0, and∑kβk=1\\sum\_\{k\}\\beta\_\{k\}=1\.

#### Joint objective \(learning matchings for tours\)\.

Training couples the spanning\-tree and matching components through the downstream tour length\. Concretely, the TSP model producesXXwhich is then decomposed throughFW Decompositionover the spanning\-tree polytope, and that yields a short list of trees with weights\{αt\}\\\{\\alpha\_\{t\}\\\}\. For each treett, we form its odd setOtO\_\{t\}, which then goes through the neural network which outputsStS\_\{t\}onEOtE\_\{O\_\{t\}\}, and we follow that by applying FW to obtain matchings\{Mt,k\}\\\{M\_\{t,k\}\\\}with weights\{βt,k\}\\\{\\beta\_\{t,k\}\\\}\. We evaluate the Christofides tour length for each\(t,k\)\(t,k\)pair and minimize the mixture objective

ℒtour=∑tαt​∑kβt,k​\(Christofides​\(Tt,Mt,k\)\),\\mathcal\{L\}\_\{\\text\{tour\}\}=\\sum\_\{t\}\\alpha\_\{t\}\\sum\_\{k\}\\beta\_\{t,k\}\\;\\\!\\left\(\\textsc\{Christofides\}\(T\_\{t\},M\_\{t,k\}\)\\right\),whereChristofides​\(Tt,Mt,k\)\\textsc\{Christofides\}\(T\_\{t\},M\_\{t,k\}\)represents the length of the TSP tour generated by the given parameters\. So the learned matchings are directly trained to be*useful for producing short tours*under our pipeline\. We also include reconstruction loss in the same way as with other problems, in order to keep X close or inside the polytope\. This training also does not require ground\-truth matchings; gradients are driven by the downstream TSP objective through the Christofides rounding\.

#### Matching network architecture\.

The learned matching utilizes a lightweight 2 layer MLP structure to perform fast inferences, it encodes nodes using Fourier positional encoding of the 2D coordinates as input features\. The matchingFW Extensionuses a small iteration budgetTmatch=8T\_\{\\text\{match\}\}=8\. The TSP neural network working on the spanning tree polytope follows the same setup as[SectionD\.1](https://arxiv.org/html/2609.25728#A4.SS1)\.

Empirically, the learned matching variant improves objective values over the greedy\-matching version \([Table10](https://arxiv.org/html/2609.25728#A4.T10)\) at the cost of higher inference time, and the runtime breakdown \([Table11](https://arxiv.org/html/2609.25728#A4.T11)\) shows that matching can dominate the additional overhead on large instances\.

### D\.6Ablation Study and Discussion on Projection

Following the experiments for Maximum Coverage, we ablate the effect of projecting the output of the neural network back into the feasible polytope\. Table[13](https://arxiv.org/html/2609.25728#A4.T13)shows that, for TSP, removing projection yields slightly better performance on both TSP\-500 and TSP\-1000\.

We attribute this to a structural property of \(2D\) Euclidean TSP\. Unlike problems with highly heterogeneous instance \(e\.g\., QAP\), Euclidean TSP is constrained by a low\-dimensional geometric embedding and metric distances, which makes extremely “hard” out\-of\-distribution instances comparatively rare—both for synthetic generators and for real\-world benchmarks that remain within the same problem family\. To support this claim empirically, we further test the model trained on TSP\-1000 on 2D Euclidean instances fromTSPLIB\[[53](https://arxiv.org/html/2609.25728#bib.bib68)\]with 51 to 1002 nodes\. As shown in[Table14](https://arxiv.org/html/2609.25728#A4.T14), our method remains robust even on this substantially broader collection\.

These results suggest that, for problems like Euclidean TSP, explicit projection is not only unnecessary but can be counterproductive, while the no\-projection design is fully viable and can even yield better performance in practice\.

Table 13:Ablation results with and without projectionTable 14:Performance on TSPLIB of the model trained on TSP\-1000\.
### D\.7Ablation Study: LP \+ Noise

To test whether the gains come from*learned*structure rather than arbitrary perturbations of the LP relaxation, we construct a randomized baseline by adding i\.i\.d\. Gaussian noise to the subtour\-LP solution:x~=xLP\+ϵ\\tilde\{x\}=x\_\{\\mathrm\{LP\}\}\+\\epsilonwithϵ∼𝒩⁡\(0,σ2\)\\epsilon\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}\)andσ=0\.2\\sigma=0\.2\. We then run the*same*FW decomposition and downstream rounding pipeline usingx~\\tilde\{x\}as the input weights\. For each instance we repeat this procedure for 10 independent noise draws and report the*best*\(shortest\) tour length, which is a favorable “best\-of\-10” setting for the noise baseline\. As shown in Table[15](https://arxiv.org/html/2609.25728#A4.T15), this randomized perturbation is unable to match our method, with the gap becoming especially pronounced on larger instances \(TSP500–TSP1000\)\. This indicates that our network learns instance\-dependent signals that improve the pipeline beyond what can be obtained by injecting random noise intoxLPx\_\{\\mathrm\{LP\}\}\.

Table 15:Ablation results on replacing the learned edge weights with random perturbations of the LP solution\. Lower is better\.

## Appendix EFurther Discussion on Inference Settings

In the field of Neural Combinatorial Optimization \(NCO\), strategies to refine solution quality*during inference*are widely adopted\. While terminology varies across domains—often referred to as Test\-Time Optimization \(TTO\) or Local Search/Sampling—these procedures share a common goal: to improve the solution beyond a single neural forward pass\. Broadly, they fall into three categories:

1. 1\.Gradient\-based Active Search:Fine\-tuning model parameters on individual test instances \(e\.g\., active search\)\.
2. 2\.Search\-based Refinement:Applying discrete search algorithms, such as beam search, MCTS, or heuristic local search \(e\.g\., 2\-opt\), to the model’s output\.
3. 3\.Parallel Sampling:Running multiple inference trajectories, such as multiple diffusion noise seeds, batched augmentations, or independent neural forward passes, to enlarge the candidate pool and select the best solution\. We distinguish this from native candidate\-set decoding, where a single neural forward pass deterministically or procedurally produces several candidates as part of the method’s prescribed decoder\.

#### Decoupling Learned Priors from Search Strategies\.

While these inference\-time techniques significantly enhance performance, they introduce a critical challenge in evaluation:they often obscure the intrinsic quality of the learned neural model\.

Recent studies have critically examined this phenomenon\.\[[6](https://arxiv.org/html/2609.25728#bib.bib65)\]demonstrated that in tree\-search\-based neural methods, the search algorithm frequently performs the majority of the "heavy lifting," rendering the contribution of the learned neural guidance marginal or, in some cases, redundant\.\[[68](https://arxiv.org/html/2609.25728#bib.bib64)\]finds that this is also true for tree\-search based heatmap recovery\. They further point out that this inconsistency is particularly concerning because the search algorithm, which is critical for determining the final solution, is not integrated during neural network training\. This confounding effect where the solver relies more on inference time compute than on the learned structural prior has also been noted in recent works such asUCOM2\[[7](https://arxiv.org/html/2609.25728#bib.bib39)\]andGeoNCO\[[23](https://arxiv.org/html/2609.25728#bib.bib11)\], which sets a clear difference between TTO and non\-TTO methods and make no comparison between them\.

#### Evaluation Protocol: The "One\-Shot" Setting\.

In this work, our primary goal is to evaluate the intrinsic constructive capability of the proposed framework\. We aim to assess how well the model aligns with the problem structure in a minimal\-latency setting, without relying on extensive post\-hoc search\. Consequently, we strictly standardize the evaluation to anon\-TTOcontext:

- •For Maximum Coverage baselines \(e\.g\.,CardNN,UCOM2,GeoNCO\), we utilize their "short" inference modes which comes without gradient\-based finetune on inference or iterative refinement steps\.
- •For TSP baselines \(e\.g\.,DIMES,DIFUSCO,COExpander\), which often rely on active search, parallelized diffusion steps or local search to refine results, we restrict them to their greedy or native candidate\-set decoding equivalents\. We exclude parallel sampling strategies that enlarge the candidate pool through repeated inference, such as running diffusion\-based solvers with multiple noise seeds\. However, when a method’s native decoder generates multiple feasible candidates from one neural network pass and selects the best one, we still regard it as one\-shot\. For example,DIMES \(RL\+S\)uses a number of samples generated from a single neural network pass, so this is treated as native candidate\-set decoding rather than parallel sampling\. Our FW decomposition follows the same principle: one neural output is decomposed into a sparse set of feasible vertices, and inference reports the best vertex among them\.

By isolating the raw performance of the underlying neural architectures, we provide a clearer benchmark of their structural learning capabilities, ensuring that improvements are attributed to better representation learning rather than increased inference computation\. This protocol also keeps the comparison consistent with the decoding budgets used by strong baselines, where candidate\-set selection is often part of the native decoder, while excluding additional inference\-time computation whose main effect is to enlarge the search pool\.

## Appendix FGPU Memory Usage

We report the GPU memory footprint of FWNCO across all training configurations\. Table[16](https://arxiv.org/html/2609.25728#A6.T16)reports peak CUDA allocated memory during training\. These measurements describe allocated GPU memory, rather than total device\-memory consumption, which can additionally include reserved but unused memory and CUDA runtime overhead\.

Table 16:Peak CUDA allocated memory during FWNCO training\.ProblemDataset /kkPeak memory \(GiB\)MCRandom500,k=10k=100\.135MCRandom500,k=50k=500\.135MCRandom1000,k=20k=200\.194MCRandom1000,k=100k=1000\.194TSPTSP500\.043TSPTSP1000\.106TSPTSP5002\.012TSPTSP10007\.708QAPQAP325\.168Overall, the framework has a relatively modest memory footprint of below 8 GiB for every configuration\. These results indicate that FWNCO does not require the latest high\-end GPUs for the evaluated configurations and suggest potential for scaling to substantially larger instances\.

## Appendix GPolytopes with Efficient \(Approximate\) LMO

#### Matroid polytope

Matroid polytopes provide a unified geometric framework for a broad class of combinatorial constraints\. Given a matroidℳ=\(V,ℐ\)\\mathcal\{M\}=\(V,\\mathcal\{I\}\)with rank functionr⁡\(⋅\)r\(\\cdot\), Edmonds’ classical formulation characterizes the associated matroid base polytope as

𝒫base=\{𝐱∈ℝ≥0n\|∑i∈Axi≤r\(A\)∀A⊆V,∑i∈Vxi=r\(V\)\}\.\\mathcal\{P\}\_\{\\text\{base\}\}=\\left\\\{\\mathbf\{x\}\\in\\mathbb\{R\}^\{n\}\_\{\\geq 0\}\\;\\middle\|\\;\\sum\_\{i\\in A\}x\_\{i\}\\leq r\(A\)\\ \\forall A\\subseteq V,\\ \\sum\_\{i\\in V\}x\_\{i\}=r\(V\)\\right\\\}\.A key property of this polytope is that it admits an efficient LMO: optimizing a linear function over𝒫base\\mathcal\{P\}\_\{\\text\{base\}\}reduces to finding a maximum\-weight base of the matroid, which can be solved in polynomial time via the greedy algorithm\. This class of polytopes subsumes several fundamental constraints as special cases:

- •Cardinality constraints \(uniform matroid\)\.The ground set isV=\{1,…,n\}V=\\\{1,\\dots,n\\\}with rank functionr⁡\(A\)=min⁡\{\|A\|,k\}r\(A\)=\\min\\\{\|A\|,k\\\}\. The corresponding base polytope is the hypersimplexΔn,k=\{𝐱∈ℝ≥0n\|∑i=1nxi=k,0≤xi≤1\}\.\\Delta\_\{n,k\}=\\left\\\{\\mathbf\{x\}\\in\\mathbb\{R\}^\{n\}\_\{\\geq 0\}\\;\\middle\|\\;\\sum\_\{i=1\}^\{n\}x\_\{i\}=k,\\;0\\leq x\_\{i\}\\leq 1\\right\\\}\.
- •Partition matroids\.LetV=⋃j=1cVjV=\\bigcup\_\{j=1\}^\{c\}V\_\{j\}be a partition of the ground set with associated budgetsk1,…,kck\_\{1\},\\dots,k\_\{c\}\. The rank function isr⁡\(A\)=∑j=1cmin⁡\{\|A∩Vj\|,kj\}r\(A\)=\\sum\_\{j=1\}^\{c\}\\min\\\{\|A\\cap V\_\{j\}\|,k\_\{j\}\\\}, and the base polytope is given by\{𝐱∈ℝ≥0n\|∑i∈Vjxi=kj​∀j\}\.\\left\\\{\\mathbf\{x\}\\in\\mathbb\{R\}^\{n\}\_\{\\geq 0\}\\;\\middle\|\\;\\sum\_\{i\\in V\_\{j\}\}x\_\{i\}=k\_\{j\}\\ \\forall j\\right\\\}\.
- •Spanning tree constraints \(graphic matroid\)\.Given an undirected graphG=\(V,E\)G=\(V,E\), the ground set consists of edges and the rank function isr⁡\(A\)=\|V\|−κ⁡\(A\)r\(A\)=\|V\|\-\\kappa\(A\), whereκ⁡\(A\)\\kappa\(A\)denotes the number of connected components in\(V,A\)\(V,A\)\. The corresponding base polytope is\{𝐱∈ℝ≥0\|E\|\|∑e∈Exe=\|V\|−1,∑e∈Axe≤\|V\|−κ\(A\)∀A⊆E\}\.\\left\\\{\\mathbf\{x\}\\in\\mathbb\{R\}^\{\|E\|\}\_\{\\geq 0\}\\;\\middle\|\\;\\sum\_\{e\\in E\}x\_\{e\}=\|V\|\-1,\\;\\sum\_\{e\\in A\}x\_\{e\}\\leq\|V\|\-\\kappa\(A\)\\ \\forall A\\subseteq E\\right\\\}\.

Our setup naturally covers all of these cases within a single algorithmic framework\. In contrast,[Karalias et al\. \[23\]](https://arxiv.org/html/2609.25728#bib.bib11)study these settings in a case\-specific manner, requiring the design of separate projection operators onto each polytope as well as bespoke decomposition procedures tailored to each constraint\. By relying solely on the availability of an efficient LMO, our results apply uniformly across all these matroidal constraints without modifying the algorithm\.

#### Birkhoff polytope

Ann×nn\\times nmatrix is called*doubly stochastic*if all of its entries are nonnegative and the sum of the entries in each row and each column equals one\.*Permutation matrices*form a special subclass of doubly stochastic matrices, characterized by having exactly one entry equal to one in each row and column and zeros elsewhere\. The set of alln×nn\\times ndoubly stochastic matrices constitutes a convex polytope known as the*Birkhoff polytope*ℬn\\mathcal\{B\}\_\{n\}\. The Birkhoff polytope lies in an\(n−1\)2\(n\-1\)^\{2\}\-dimensional affine subspace ofn2n^\{2\}\-dimensional Euclidean space defined by2​n−12n\-1independent linear constraints specifying that the row and column sums all equal 1\. The extreme points of Birkhoff polytopeℬn\\mathcal\{B\}\_\{n\}are exactly the permutation matrices\. Given a residual matrixR∈ℝn×nR\\in\\mathbb\{R\}^\{n\\times n\}, the LMO overℬn\\mathcal\{B\}\_\{n\}solves

P⋆∈arg​maxP∈ext⁡\(ℬn\)⁡⟨R,P⟩=∑i∑jR⁡\(i,j\)​P​\(i,j\)P^\{\\star\}\\;\\in\\;\\operatorname\*\{arg\\,max\}\_\{P\\in\\operatorname\{ext\}\(\\mathcal\{B\}\_\{n\}\)\}\\langle R,P\\rangle=\\sum\_\{i\}\\sum\_\{j\}R\(i,j\)P\(i,j\)This problem is equivalent to finding a maximum\-weight perfect matching in a complete bipartite graphG=\(U,V,E\)G=\(U,V,E\), whereUUandVVindex the rows and columns ofRR, respectively, and each edge\(i,j\)∈E\(i,j\)\\in Ehas weightR⁡\(i,j\)R\(i,j\)\. Hence the linear maximization problem over Birkhoff polytope can be solved using the Hungarian algorithm in timeO⁡\(n3\)O\(n^\{3\}\)\. Hence, our approach subsumes the setting considered in\[[50](https://arxiv.org/html/2609.25728#bib.bib21)\]without requiring any projection onto the Birkhoff polytope\. In contrast, the authors of\[[50](https://arxiv.org/html/2609.25728#bib.bib21)\]develop a specialized decomposition algorithm and introduce an additional loss term to penalize neural network outputs that lie outside the polytope\.

#### DMO: greedy bipartite matching\.

Note that a simple greedy algorithm produces a 0\.5\-approximation to maximum weight matching for any \(positively\) weighted graph\. This still outputs a rather good matching while being significantly faster than running an exact Hungarian algorithm oracle at every iteration\. Empirically, we observe that the greedy DMO does not sacrifice quality\.

#### Matching polytope

Building on the discussion of the Birkhoff polytope, our approach also applies to matching and perfect matching polytopes\.