Combinatorial Optimization with Pointer Networks: Theory, Benchmarks, and Limitations
Intro
This article is the opening of an article series that, hopefully, should guide us and help me personally dive into the world of Machine Learning for combinatorial tasks.
Let me start by introducing the type of problems we are about to discuss further.
Combinatorial Optimization Problems
COPs are problems where the goal is to find an optimal solution from a finite set of feasible solutions.
Formal Math Definition (Ausiello et al., Erwig & Kumar)
Formally, a combinatorial optimization problem can be defined as a quadruple:
where:
- I: The set of all possible problem instances.
- f(x): For a given instance x∈I, this function maps to the finite set of feasible solutions (choices that satisfy all rules and constraints).
- m(x,y): A measure function (or objective function) that assigns a numerical value/cost to a feasible solution yy.
- g: The optimization goal, which is either to minimize or maximize the value of m(x,y).
The ultimate goal is to find an optimal solution y∈f(x) whose score achieves the optimum:
A Quick Concrete Example: The Travelling Salesman Problem
To see how this maps to reality, consider the Travelling Salesman Problem (TSP):
- Instance x∈I: A specific set of n cities with their 2D coordinates.
- Feasible solutions f(x): All valid permutations of cities that form a single closed loop visiting every city exactly once and returning to the start.
- Solution y: One specific tour (e.g., City 1 → City 4 → City 2 → City 3 → City 1).
- Measure m(x,y): The total tour length (e.g., sum of Euclidean distances between consecutive cities).
- Goal g: min (find the tour with the shortest total distance).
In practical solver formulations (such as Integer Linear Programming), a solution y is represented through a set of decision variables (in the example from the above, binary variables xij ∈ {0,1} indicating whether the path travels directly from city i to city j). The problem’s constraints (e.g., each city must have exactly one entry and one exit, without isolated sub-loops) filter out invalid configurations and define the feasible set f(x). The goal is to find a decision variable assignment that satisfies all constraints while minimizing m(x,y)m(x,y).
SOTA Production Solutions
Let’s assume you have an online shop (or delivery service) and you want to coordinate vehicles to deliver goods to customers. Or you have a warehouse and a set of orders that look like "go somewhere, pick something up, and put it somewhere else." Or you have a PCB where it is required to drill thousands of tiny holes (TSP for PCB). You want to solve these problems optimally according to your needs, e.g., minimize the travel time of the vehicles or minimize and evenly distribute the battery consumption of the robots in the warehouse (VRP problem).
In practice, how does industry tackle these problems today? Existing production solutions generally fall into two broad categories:
1. Exact Solvers (Branch-and-Bound / MIP)
Solvers like Concorde (the gold standard for TSP) or mathematical integer programming engines like Gurobi and CPLEX search systematically through the solution space using techniques like branch-and-cut.
- The upside: They provide mathematical guarantees of optimality.
- The bottleneck: Because these problems are NP-hard, the worst-case time complexity grows exponentially (O(2ⁿ)). While Concorde can find provably optimal solutions for dozens or hundreds of cities, it quickly becomes computationally intractable when instances grow or when latency budgets are tight.
2. Heuristics, Metaheuristics & Production Solvers
In many scenarios, the time required to find a solution is so critical that an exact solution is unfeasible and a high-quality suboptimal solution is more than enough. This is the domain of Approximation Algorithms (Wikipedia), local search, and industrial optimization suites:
- Google OR-Tools (OR-Tools): Google’s widely used open-source suite for routing and constraint programming. It pairs construction heuristics (e.g., Cheapest Arc) with local search metaheuristics (such as Guided Local Search, Tabu Search, or Simulated Annealing) to navigate complex vehicle routing and scheduling constraints.
- Hexaly (formerly LocalSolver, Hexaly): A state-of-the-art commercial solver that combines exact methods with high-performance local search heuristics, designed specifically to tackle large-scale vehicle routing and packing problems with tight turnaround times.
- Specialized Heuristics (e.g., LKH-3): Iterative heuristics like Lin-Kernighan (k-opt edge swapping) routinely produce solutions within 0.5% of optimum even for thousands of nodes.
These production solvers are powerful, mature, and handle complex real-world operational scenarios. However, they share a fundamental characteristic, which sometimes is a drawback: they solve every new instance from scratch via iterative search. For each incoming query, the solver must perform thousands or millions of neighborhood transitions, which can take up several minutes (or even hours) per problem.
Why Machine Learning in Combinatorial Optimization?
In modern latency-critical environments—such as real-time ride dispatching, dynamic robot pathing, or on-the-fly request rerouting—waiting several seconds or minutes for a solver loop to converge is often a deal-breaker.
This is where Machine Learning offers a radically different paradigm:
- Near-Optimal is Good Enough: In the real world, physical execution is noisy (traffic jams, loading delays, weather). A solution within 2%–4% of the global optimum delivered instantly is far more valuable than a provably optimal tour that arrives too late.
- Amortized Computation: Instead of spending heavy compute on iterative search at runtime, we front-load the computational cost into an offline training stage. Over millions of synthetic problem instances, a neural network learns the general spatial and relational heuristics of what makes a tour efficient.
- Sub-Millisecond Inference & Massive Scalability: At inference time, finding a solution requires only a single forward pass through the network (or greedy autoregressive decoding). On modern GPUs, this takes fractions of a millisecond per instance and can easily evaluate batches of hundreds of instances in parallel.
So the question is, how do we design a neural network capable of taking a variable-sized set of city coordinates and directly outputting an optimal permutation of them?
Pointers Network
Here I will mostly try to explain in the simple terms this Pointer Network article by Oriol Vinyals, Meire Fortunato, and Navdeep Jaitly with small additional comments and adjustments from my end.
To make it easier to understand, we will start with the TSP problem. The general objective we would have in this case would be to generate the sequence of the cities the salesman has to visit in order to travel the least possible distance.
Allow me to split the task into iterations. Let’s say we have 7 cities in total and the salesman is currently standing at the 3rd one, assuming the first two cities were selected in the most optimal way. Now, suppose we could somehow pack the overall geographical information of all the cities (for simplicity, seven 2D vectors [x, y]) along with the trace of the previously visited cities to decide where to move next. How do we do that? We can estimate the likelihood of each candidate city j being the most optimal next choice.
Here, aᵗ is a likelihood. It is estimated as the output of the softmax function, which is the standard way to convert raw estimated scores into a valid probability distribution:
where:
uⱼᵗ — the scalar score (logit) of city j at step t (for already visited cities, we mask their scores to −∞ so their probability becomes 0).
exp(uⱼᵗ) (or e^(uⱼᵗ)) — the exponential of the score for city j (makes all values strictly positive and sharply amplifies the leader’s margin).
∑ exp(uₖᵗ) — the sum of exponentials of the scores across all n cities (the normalizing denominator that ensures all probabilities sum to one: ∑ aⱼᵗ = 1).
aⱼᵗ ∈ (0, 1) — the final probability of selecting city j at step t. Afterwards, we greedily pick the city with the highest probability, assuming that the city it points to is our optimal next step:
In our example, this points to City 5 with a probability of 0.68.
Now, let’s figure out how the logit uⱼᵗ is calculated.
Get Anton Elphimov’s stories in your inbox
Join Medium for free to get updates from this writer.
How would you find the probability that the tree in front of you is the tallest in the neighborhood? The natural answer is to estimate the height of that tree and compare it against the height estimates of all the other trees. That is essentially how things work in machine learning as well, with one additional twist: we also need to learn the reference ground level, since it is not given in advance. The only difference is that our “ground” is multidimensional, making it a hyperplane rather than a flat surface.
Assume we have extracted a combined feature vector sⱼᵗ that holds both the information about candidate city j and the trace of previously visited cities up to step t (imagine this point as the top of the tree). We also have a learnable vector v, which serves as the normal vector to the ground hyperplane in our feature space. To estimate the score (the height of the tree), we simply compute the projection of the feature vector sⱼᵗ onto the normal vector v (I encourage you to check how the distance between a point and a plane is measured if you feel rusty):
P.S. Since we are only interested in relative comparison values, normalizing by the magnitude ‖v‖ can be safely omitted here.
The question now is, where do we get this feature vector sⱼᵗ from? As mentioned earlier, it combines information about the candidate city we are currently checking (let’s call this vector eⱼ) and information about the cities visited so far (let’s call it dₜ). I don’t know about you, but when I need to combine two things, my first instinct is to take a weighted sum:
W₁eⱼ + W₂dₜ.
Can we use it directly? Not quite. If we left it like that, the second term W₂dₜ would remain identical across all candidate cities j, since the traversed path is already fixed at step t. That means our score would split into two independent parts:
And when plugged into the softmax, that constant C simply cancels out
Simply put, in a purely linear setup, W₂dₜ would play no role at all! Yet intuition tells us otherwise: if you are currently in Berlin, it clearly makes much more sense to visit Paris first and then fly to Washington than the other way around, right?
So, to truly fuse both parts together and prevent them from separating, we wrap the sum in a non-linear activation function — specifically, the hyperbolic tangent (tanh). This gives us our feature vector sⱼᵗ = tanh(W₁eⱼ + W₂dₜ) and the complete attention mechanism:
To make the picture complete, we need to understand how to obtain eⱼ and dₜ. Let’s focus on eⱼ, which represents the feature vector of a city. The task is to convert the geographical coordinates of the cities into a collection (matrix) of feature vectors in a fixed-dimensional space. At the time the Pointer Network paper was published, the dominant architecture for sequence modeling was the LSTM (Long Short-Term Memory).
Before Transformers took over, LSTMs were the primary architecture used to encode sequential data (such as words in a sentence or coordinates in a tour) into rich feature vectors.
You can think of an LSTM as a recurrent cell rolled out over time. At each step t, the cell takes two primary inputs: the current input vector (xₜ, e.g., city coordinates or a word embedding) and the previous hidden state (hₜ₋₁). Rather than simply adding them together, the LSTM maintains a long-term memory highway called the cell state (cₜ) alongside the working hidden state (hₜ). It uses internal gating mechanisms — Forget, Input, and Output gates — to control how much information to erase, add, and emit at each step
This gating mechanism allows the network to capture dependencies across long sequences. In the Pointer Network, all input cities are processed by this encoder to build the sequence of representations e₁, …, eₙ. In the same way, the feature vector for the cities visited so far is constructed using the LSTM decoder.
And now we have everything in place and ready for the final scheme.
Tests & Experimental Results
To see how the Pointer Network performs in practice, I trained and evaluated the architecture on the Travelling Salesman Problem using the combinatorics-ml repository across different scales: fixed small size (N = 20), medium size (N = 100), and a multi-size variable scale regime (N ∈ [100, 200]).
Experimental Setup
- Environment: 2D Euclidean TSP with coordinates sampled uniformly from the unit square 𝒰(0, 1)².
- Datasets: 128,000 training instances and 1,000 validation instances per size.
- Architecture: 2D coordinate embedding (dim 128), bidirectional LSTM encoder (hidden dim 128), and LSTM decoder with additive pointer attention and dynamic masking of visited cities.
- Training: REINFORCE policy gradient with an Exponential Moving Average (EMA, β = 0.9) baseline, trained for 10 epochs on a Tesla T4 GPU.
- Baseline Solver: Google OR-Tools using
PATH_CHEAPEST_ARCfollowed byGUIDED_LOCAL_SEARCHwith a 10-second per-instance time limit.
The relative quality of the network is quantified using the optimality gap against OR-Tools:
1. Scaling from TSP-20 to TSP-100
Evaluating the trained checkpoints with greedy decoding against Google OR-Tools
TSP-20 Tour Examples
On small instances (N = 20), the network learns near-perfect geometric intuition. It connects nearby nodes in natural convex perimeters without edge crossings:
TSP-100 Tour Examples
Scaling the same architecture directly to 100 cities (N = 100) still yields coherent macro-tours in under half a millisecond per instance:
While the network avoids gross self-intersections and traces clean global clusters, you can spot occasional sub-optimal crossings on denser clusters that an iterative local search would untangle.
2. Multi-Size & Zero-Shot Generalisation (N = 100 … 200)
Because the attention pointer mechanism is independent of sequence length, a single Pointer Network can be evaluated on arbitrary problem sizes. To test variable-size scaling and zero-shot generalisation, a model was trained on a mixture of sizes (N ∈ {100, 120, 140, 160}) and then tested on both seen sizes and unseen holdout sizes (N ∈ {150, 180, 200}):
The plots below illustrate the validation trajectories, tour length scaling, and the optimality gap across problem sizes:
And here is the model generating zero-shot solutions on seen (N = 160) and unseen problem sizes (N = 150, 180, 200):
Key Takeaways & Discussion
1. Solution Quality Decays with Problem Size
Looking across the spectrum — from N = 20 up to N = 200 — a clear pattern emerges:
- At TSP-20, the Pointer Network achieves an impressive +4.06% gap relative to Guided Local Search.
- At TSP-100, the gap expands to +7.55% ~ +8.57%.
- At TSP-200, the gap climbs to +10.88%.
The larger the problem instance, the greater the combinatorial complexity, and the harder it becomes for a purely greedy autoregressive pointer to avoid suboptimal local commitments early in the tour.
2. The Fundamental Flaw: Is LSTM Really the Right Encoder?
The Pointer Network was a groundbreaking milestone because it demonstrated that neural networks could output combinatorial permutations. However, from a modern architectural standpoint, using an LSTM as the encoder reveals a severe mismatch with the problem:
- A TSP instance is a Set, not a Sequence:
A collection of 2D cities {x₁, x₂, …, xₙ} has no inherent temporal or sequential order. Whether city A is passed to the network first or fifth should not change its mathematical identity or its spatial relationships with its neighbors. - Order Sensitivity & Recency Bias:
Because an LSTM processes cities recurrently (x₁ → x₂ → … → xₙ), it imposes an artificial sequential bias. The final hidden state hₙ has compressed x₁ through n recurrent gates, whereas xₙ is fresh in memory. Permuting the arbitrary input order of cities leads to completely different embeddings! - The Transformer / Self-Attention Alternative:
This strongly suggests that a Transformer encoder (Multi-Head Self-Attention) is fundamentally better suited for combinatorial problems. Multi-head self-attention operates natively on sets: it is permutation-equivariant, gives every node an 𝒪(1) direct connection to every other node, and models spatial pairwise geometry without recurrent decay.
In the upcoming parts of this series, we will leave recurrent architectures behind and explore how Transformer-based Attention Models (AM) and Graph Neural Networks (GNNs) resolve these encoder bottlenecks to achieve tighter optimality gaps on larger combinatorial tasks.
All code, data generation scripts, and baseline evaluation notebooks are open-sourced in the combinatorics-ml repository.
