Heuristic — constructive

▶ Open interactive explainer →

Fourier-basis Constructive Solver
Alias fourier
Type Heuristic — constructive
Complexity O(K_max · epochs · (n + M) log M) per run

Description

Encodes a TSP tour as a closed curve in the complex plane and optimises the Fourier coefficients of that curve with gradient descent. Decoding is a pure argsort: no penalty, no repair step, and no possibility of producing an invalid tour.

Curve representation

The tour is the closed curve

γ(s) = Σ_{k=-K}^{K} c_k · exp(2πi · k · s),   s ∈ [0, 1)

sampled at M evenly-spaced points s_j = j/M. The 2K+1 complex coefficients c_k are the only free variables; gradient descent moves them so the curve passes near each city.

Energy

E(c) = Σ_i  min_j |city_i − γ(s_j)|²   +   λ Σ_k (2πk)² |c_k|²
         attraction                              tension
  • Attraction: each city pulls its nearest curve sample toward it.
  • Tension: a smoothness regulariser that weights high-frequency modes proportionally to their squared Fourier frequency, preventing jagged curves that visit no city cleanly.

Coarse-to-fine optimisation loop

initialise c[0] = centroid, c[1] = radius/2, all others = 0
λ ← opts.lambda

for k_active = 1 … K_max:
    basis[k][j] = exp(2πi · ks[k] · j/M)   // pre-computed; constant within stage

    repeat opts.epochs times:
        γ = eval_curve(c, ks, M)
        grad = attraction_gradient + tension_gradient
        for k where |ks[k]| ≤ k_active:
            c[k] -= (lr / n) * grad[k]

    λ *= lambda_decay

Unlocking modes one stage at a time lets the optimiser set the overall loop shape (low modes) before refining local detail (high modes), avoiding the saddle points that occur when all modes compete simultaneously.

Decode (always valid)

s_i = argmin_j |city_i − γ(s_j)|     // nearest curve sample per city
tour = argsort(s_i)                    // sort cities by their curve position

The argsort gives array positions in cities[], which are then mapped to their .id fields — the same pattern used by Christofides. This guarantees a valid Hamiltonian tour regardless of convergence quality.

procedure FourierSolver(cities):
    coeffs ← fit_fourier_basis_to_cities(cities)
    curve ← reconstruct_closed_curve(coeffs)
    tour ← argsort_cities_by_curve_parameter(cities, curve)
    tour ← two_opt_polish(tour)
    return tour

Options

Field CLI flag Default Range Description
k_max --k-max 4 ≥ 1 Maximum Fourier mode (number of frequency stages)
m --m 200 ≥ 2 Curve sampling resolution (points on γ)
lambda — 0.05 > 0 Initial tension weight
lambda_decay — 0.5 (0, 1) Tension multiplier applied at each stage
lr — 0.05 > 0 Gradient descent learning rate
epochs --epochs 400 ≥ 1 Gradient steps per k_active stage

epochs follows the same vocabulary as all other solvers in this codebase, with one exception: unlike HeuristicOptions.epochs elsewhere, epochs=0 is not a “run forever” sentinel here — it’s rejected by validation, since it counts gradient steps per k_active stage, not outer solver iterations.

Only k_max, m, and epochs have CLI flags; lambda, lambda_decay, and lr are reachable only via the REST API’s configs.fourier or a [fourier]/[stage.fourier] TOML table (field names match the table above).

Tuning for larger instances

k_max (coefficient count) is the dominant quality lever, not m (curve resolution). On a280 (280 cities), scaling k_max alone took the standalone gap from +103.4% at the shipped default (k_max=4) down to +24.6% at k_max=32 (further down to +7.6% after piping into a 2-opt polish pass, which is the more relevant number since Fourier is typically used as a warm-start). Scaling m alone, without more coefficients, made quality worse.

Config Standalone gap +2-opt gap Wall time
k_max=4, m=200 (default) +103.4% — 0.19s
k_max=32, m=200 +24.6% +7.6% 4.25s
k_max=32, m=1120 (4n) +15.9% +15.8% 19.2s

(a280, 280 cities; wall times measured after the KD-tree nearest-sample optimization that shipped in #370.)

Cost scales as k_max × epochs gradient steps total, with no upper bound enforced by validate() — there’s no instance-size-aware guard, so an oversized --k-max on a large instance will simply run for a long time (or, at extreme values, effectively hang) rather than error out. Start by scaling k_max alone before touching m; m scaling is comparatively expensive and, per the table above, doesn’t reliably help quality on its own.

Usage

# standalone
teeline solve fourier -i ./data/tsplib/berlin52.tsp

# as warm-start for 2-opt (recommended)
teeline pipeline --steps=fourier,2opt -i ./data/tsplib/berlin52.tsp

# as warm-start for LK
teeline pipeline --steps=fourier,lk -i ./data/tsplib/berlin52.tsp

# tuning k_max for a larger instance
teeline solve fourier -i ./data/tsplib/a280.tsp --k-max 32

Per-stage TOML config (via pipeline --config):

[[stage]]
solver = "fourier"

[stage.fourier]
k_max = 32
m     = 200

Notes

  • init_tour is ignored: as a purely constructive solver, Fourier does not accept or benefit from a seed tour; the --init-tour pipeline option has no effect on this solver.
  • f64 internally: all gradient computation uses f64 for numerical stability; coordinates are converted from f32 at input and the final tour is Vec<usize> city IDs.
  • Nearest-sample lookup uses a KD-tree: finding each city’s nearest point on the curve (used in both the attraction gradient and the final decode) queries a KD-tree built fresh from the current curve samples each step, rather than a linear scan. The KD-tree module stores coordinates as f32, so the which-sample-is-nearest decision is made at f32 precision; the gradient and decode math around the returned index still run in f64. The query point is tagged with a sentinel id (usize::MAX) that can never collide with a real curve-sample id, since the KD-tree’s self-exclusion guard (designed for querying a point against the same id space it came from, as in nearest_neighbor.rs) would otherwise make curve sample 0 permanently unselectable whenever the query’s id defaulted to 0.
  • WASM-compatible: uses num-complex = "0.4" (no_std-compatible, no libm dependency).

Relationship to the Elastic Net

This algorithm is a Fourier-parameterised variant of the Elastic Net (Durbin & Willshaw 1987). Both share the same two-term energy (city attraction + curve smoothness/tension) and optimise via gradient descent. The key differences:

Elastic Net This implementation
Curve representation M explicit node positions 2K+1 Fourier coefficients
Tension term Sum of squared edge lengths λ(2πk)² |c_k|² (diagonal in coefficient space)
Mode schedule Simultaneous + temperature annealing Coarse-to-fine frequency unlocking
Tour decode Explicit node ordering Argsort of nearest-sample parameter s_i

References

  • Durbin, R. & Willshaw, D. (1987) — “An analogue approach to the travelling salesman problem using an elastic net method”, Nature 326, 689–691. DOI: 10.1038/326689a0 (original Elastic Net; this solver is a Fourier-parameterised variant)