Exact solvers¶

Three solvers that certify the optimum and set is_optimal_. They are the yardstick of the library: every heuristic is tested against BruteForce on small instances, and MILP proves the published optima of the small Waterloo instances.

Solver	Objective it certifies	Cap (`max_nodes`)	Typical time
`BruteForce`	plain tour and multi-trip, under `split="greedy"` or `"optimal"`	11	a fraction of a second at n = 11
`HeldKarp`	plain tour, symmetric or asymmetric	20	under a second at n = 20 (~90 MB of tables)
`MILP`	plain tour, symmetric or asymmetric	300	wi29 and dj38 in well under a second; qa194 proved in about a minute

HeldKarp and MILP raise ValueError when max_time_work is given: an "exact" result that is not optimal for the objective would be a trap, so a certified multi-trip optimum comes only from BruteForce (n ≤ 11) and larger multi-trip instances go to the budget-aware heuristics.

Under split="greedy", BruteForce is exact over greedy-decoded giant tours (a partition that closes a trip while the next node would still fit is not representable); under split="optimal" it is exact for the distance-constrained multi-trip problem.

>>> from skroute import BruteForce, HeldKarp, MILP
>>> from skroute.datasets import load_tsp
>>> wi = load_tsp("wi29")
>>> est = MILP().fit(wi.distance_matrix(), labels=wi.labels)
>>> int(est.cost_), est.is_optimal_, est.gap_
(27603, True, 0.0)

BruteForce ¶

BruteForce(max_nodes: int = 11)

Bases: BaseRouter

Exact solver by exhaustive enumeration of the (n - 1)! giant tours.

Every permutation of the non-depot nodes is priced with the problem's own decoder, so the optimum is certified for the plain tour and for the multi-trip objective under the chosen split rule (fit(..., split="greedy" | "optimal")). It is the reference every other solver is tested against on small instances.

Parameters:

Name	Type	Description	Default
`max_nodes`	`int`	Hard cap on the number of nodes (`fit` raises `ValueError` above it). The work grows as `(n - 1)!`: `10!` at n = 11 is 3.6 million tours — still a fraction of a second (a few hundredths on a laptop); n = 12 is eleven times more, n = 13 a hundred and thirty. Raise it only if you accept the time.	`11`

Attributes:

Name	Type	Description
`is_optimal_`	`bool`	Always `True`: the search is exhaustive.

Notes

The permutations of positions 1..n-1 are visited in lexicographic order with Knuth's next permutation (TAOCP 7.2.1.2, Algorithm L), inside a nogil loop, and the first strictly better tour is kept, so a tie resolves to the lexicographically first optimum — the same tour itertools.permutations would find first. On a symmetric matrix without a budget the reversed tour costs the same, so permutations with tour[1] > tour[n-1] are skipped and the kept orientation is the lexicographically smaller one; under a budget every orientation is priced because the split depends on direction.

Under split="greedy" the result is exact over greedy-decoded giant tours (a partition that closes a trip while the next node would still fit is not representable); under split="optimal" it is exact for the distance-constrained multi-trip problem.

Complexity O((n - 1)! · n) time (roughly halved when symmetric and unbudgeted: the reversal of a pair is also priced when the kept orientation is within 1e-9 of the incumbent), O(n) memory.

Supports: symmetric and asymmetric matrices, multi-trip objective (both split rules); deterministic.

References

.. [1] D. E. Knuth, The Art of Computer Programming, Vol. 4A, §7.2.1.2 "Generating all permutations", Algorithm L.

Examples:

A plain TSP from numpy; the optimum is certified:

>>> import numpy as np
>>> from skroute import BruteForce
>>> C = np.array([[0, 5, 9, 10], [5, 0, 4, 8], [9, 4, 0, 3], [10, 8, 3, 0]], dtype=float)
>>> bf = BruteForce().fit(C)
>>> bf.tour_.tolist(), bf.cost_, bf.is_optimal_
([0, 1, 2, 3], 22.0, True)

The multi-trip objective from a dict-of-dicts (the depot is the first key): a 4-hour working day and a 3.0 charge per extra trip.

>>> cost = {
...     1: {1: 0, 2: 5, 3: 9, 4: 10},
...     2: {1: 5, 2: 0, 3: 4, 4: 8},
...     3: {1: 9, 2: 4, 3: 0, 4: 3},
...     4: {1: 10, 2: 8, 3: 3, 4: 0},
... }
>>> hours = {
...     1: {1: 0, 2: 1, 3: 2, 4: 2},
...     2: {1: 1, 2: 0, 3: 1, 4: 2},
...     3: {1: 2, 2: 1, 3: 0, 4: 1},
...     4: {1: 2, 2: 2, 3: 1, 4: 0},
... }
>>> bf = BruteForce().fit(cost, time_matrix=hours, max_time_work=4.0, extra_cost=3.0)
>>> bf.route_.tolist(), bf.cost_, bf.n_trips_
([1, 2, 3, 1, 4, 1], 41.0, 2)
>>> bf.trip_times_.tolist()
[4.0, 4.0]

HeldKarp ¶

HeldKarp(max_nodes: int = 20)

Bases: BaseRouter

Exact TSP solver by the Held-Karp dynamic programme over subsets of nodes.

dp[S][j] is the cheapest path from the depot through exactly the node set S ending at j; the optimum is the best dp[all][j] + C[j, depot]. Every arc is read directionally, so asymmetric matrices are solved exactly. The plain tour is the only objective it certifies: under a budget it raises (D6) — use BruteForce for a certified multi-trip optimum, or a heuristic solver.

Parameters:

Name	Type	Description	Default
`max_nodes`	`int`	Hard cap on the number of nodes (`fit` raises `ValueError` above it). Time is O(2^(n-1) · n²) and memory `2^(n-1) · (n-1)` doubles plus as many `int8` parent bytes: about 90 MB (80 MB of doubles) and well under a second at n = 20; ~830 MB at n = 23 (738 MB of doubles plus a 92 MB parent table). Raise it only if you accept the memory; `MILP` solves these sizes instantly in any case.	`20`

Attributes:

Name	Type	Description
`is_optimal_`	`bool`	Always `True`: the programme is exhaustive over subsets.

Notes

Push formulation in a nogil loop over the subsets in increasing numeric order (every proper subset of S is numerically smaller, so dp[S] is final when expanded); malloced tables freed before returning; the parent table is int8. Ties are broken by node index at each expansion, not lexicographically over whole tours as BruteForce does — equal-cost optima may come out in a different order.

Complexity O(2^(n-1) · n²) time, O(2^(n-1) · n) memory.

Supports: symmetric and asymmetric matrices; plain TSP only (raises under a budget); deterministic.

References

.. [1] M. Held and R. M. Karp, "A dynamic programming approach to sequencing problems", Journal of the Society for Industrial and Applied Mathematics 10 (1962) 196-210. .. [2] R. Bellman, "Dynamic programming treatment of the travelling salesman problem", Journal of the ACM 9 (1962) 61-63.

Examples:

>>> import numpy as np
>>> from skroute import HeldKarp
>>> C = np.array([[0, 5, 9, 10], [5, 0, 4, 8], [9, 4, 0, 3], [10, 8, 3, 0]], dtype=float)
>>> hk = HeldKarp().fit(C)
>>> hk.cost_, hk.is_optimal_
(22.0, True)
>>> int(hk.route_[0]) == int(hk.route_[-1]) == 0 and sorted(hk.tour_.tolist()) == [0, 1, 2, 3]
True

An asymmetric matrix is read directionally (the reversed tour costs 25 here):

>>> A = np.array([[0, 1, 9, 9], [9, 0, 1, 9], [9, 9, 0, 1], [1, 9, 9, 0]], dtype=float)
>>> HeldKarp().fit(A).tour_.tolist()
[0, 1, 2, 3]

MILP ¶

MILP(time_limit: float | None = 60.0, max_nodes: int = 300, mip_rel_gap: float = 0.0)

Bases: BaseRouter

Exact TSP solver: the Dantzig-Fulkerson-Johnson integer programme with lazy subtour cuts.

The tour is modelled with binary edge variables (symmetric matrix, n(n-1)/2 of them, degree-2 equalities) or arc variables (asymmetric, n(n-1), in/out-degree-1 equalities) and solved with HiGHS through scipy.optimize.milp. Subtour-elimination constraints are added lazily: after each solve the support of the integral solution is split into connected components and one cut sum(x[e] for e inside S) <= |S| - 1 per component is appended, until the solution is a single Hamiltonian cycle — a proven optimum. The plain tour is the only objective it certifies: under a budget it raises (D6).

Parameters:

Name	Type	Description	Default
`time_limit`	`float or None`	Wall-clock budget in seconds for the whole cut loop; every solve receives the time that remains. `None` runs until proven optimal. When the budget runs out the fit still returns a valid tour (see `is_optimal_`). HiGHS checks its limit only between internal phases, so a single solve can overrun the time it was handed by the length of one uninterruptible phase — its presolve above all: on the 89 700-variable asymmetric model at n = 300 a 1 s budget took 1.8 s. `fit_time_` may therefore exceed `time_limit` by roughly the duration of one solve of the model at hand.	`60.0`
`max_nodes`	`int`	Hard cap on the number of nodes (`fit` raises `ValueError` above it). Realistic sizes within a minute: ~200 symmetric nodes, ~60 asymmetric.	`300`
`mip_rel_gap`	`float`	Relative MIP gap HiGHS accepts on each solve (`0.0` = solve every relaxation to proven optimality). A positive value speeds the loop up but the final tour is then only guaranteed within that gap: `lower_bound_` comes from the dual bound and `is_optimal_` is `True` only when `gap_` is zero.	`0.0`

Attributes:

Name	Type	Description
`is_optimal_`	`bool`	`True` when the last solve returned an integral single-component solution that HiGHS reported optimal (then `lower_bound_ == cost_` and `gap_ == 0.0`); `False` on time-out. HiGHS proves optimality to absolute tolerances (`mip_feasibility_tolerance` and `mip_abs_gap`, both 1e-6, which `scipy.optimize.milp` does not expose), so the programme is solved on a normalised objective — the costs multiplied by the power of two that brings the largest one into `[8192, 16384)` — and the certificate holds to about 1e-10 times the largest cost, below the library's 1e-9 relative tolerance, whatever the units of the input (see Notes).
`lower_bound_`	`float`	The largest valid lower bound on the optimum seen: the trivial assignment bound to start with (`sum_i min_j C[i, j]`; half the two smallest entries per row when symmetric), the objective of every solve HiGHS finished (`status == 0` — a DFJ relaxation solved to optimality is a valid bound), and the dual bound `mip_dual_bound` of a time-limited solve when HiGHS reports one. The primal value of a time-limited incumbent is an upper bound and is never used. Clipped to `cost_` (a bound above a known tour is float noise).
`gap_`	`float`	`max(0.0, (cost_ - lower_bound_) / cost_)`; `0.0` when optimal.
`n_solves_`	`int`	Number of MILP solves run.
`n_cuts_`	`int`	Number of subtour-elimination constraints added.

Notes

On time-out the returned tour is the best single-component integral incumbent seen or, failing that, the core's nearest-neighbour tour polished by 2-opt (two_opt_descent when symmetric; the full-evaluation local_search_generic with the 2-opt move when asymmetric, because O(1) reversal deltas are exact only on symmetric matrices).

The objective handed to HiGHS is C multiplied by a power of two that brings its largest entry into [2**13, 2**14); every objective value and bound is divided by the same factor on the way back. The scaling is exact (a power of two rounds nothing), so the programme, its optimum and the order of tours are exactly those of the unscaled problem, but HiGHS's absolute tolerances — 1e-6 on integrality, on the primal-dual gap and on the cutoff that prunes a node "not better than the incumbent" — then bite at about 1e-10 of the largest cost whatever the units. Without it, on costs of order 1e-5 or below HiGHS prunes the true optimum as not better than its incumbent by 1e-6 and returns the runner-up with an "optimal" status, i.e. a false certificate.

DFJ is preferred to the Miller-Tucker-Zemlin formulation because MTZ's relaxation is weak: it takes hours around n = 200 where DFJ proved qa194 (9352) in about 40 s. Components are found with scipy.sparse.csgraph.connected_components (weak connectivity); HiGHS runs single-threaded and deterministically, so two fits with the same input give the same tour unless the time limit intervenes.

Progress (status, objective, number of components and bound of every solve) is logged at DEBUG level on the skroute logger.

Complexity: NP-hard in the worst case; in practice a few dozen solves of an LP-sized programme. Memory O(n²) for the variables and constraints.

Callback events (D30, D31): "start" has no tour; every cut round emits one "iteration" (its index is the solve number) whose extra carries edges (the support of the solution as (label, label) pairs), edge_weights (the solver's value of each support variable, parallel to edges and clipped to [0, 1] — 1.0 up to the solver's integrality tolerance), n_components, lower_bound, objective (the solve's objective in the units of C) and n_cuts (constraints added so far); tour is None until the support is a single Hamiltonian cycle. A callback returning True ends the cut loop like a time-out (is_optimal_ is False unless that very round had already proven the tour).

Supports: symmetric and asymmetric matrices; plain TSP only (raises under a budget); deterministic.

References

.. [1] G. Dantzig, R. Fulkerson and S. Johnson, "Solution of a large-scale traveling-salesman problem", Journal of the Operations Research Society of America 2 (1954) 393-410. .. [2] Q. Huangfu and J. A. J. Hall, "Parallelizing the dual revised simplex method", Mathematical Programming Computation 10 (2018) 119-142 (HiGHS).

Examples:

>>> import numpy as np
>>> from skroute import MILP
>>> C = np.array([[0, 5, 9, 10], [5, 0, 4, 8], [9, 4, 0, 3], [10, 8, 3, 0]], dtype=float)
>>> est = MILP().fit(C)
>>> est.cost_, est.is_optimal_, est.gap_, est.lower_bound_ == est.cost_
(22.0, True, 0.0, True)
>>> est.n_solves_ >= 1 and est.n_cuts_ >= 0
True

Western Sahara (29 cities) is proven optimal in a few seconds:

>>> from skroute.datasets import load_tsp
>>> wi = load_tsp("wi29")
>>> est = MILP().fit(wi.distance_matrix(), labels=wi.labels)
>>> int(est.cost_) == wi.optimal_tour_length == 27603, est.is_optimal_
(True, True)