Self-organising map¶
skroute.metaheuristics.SOM is the neural approach to the TSP: a one-dimensional
Kohonen ring of neurons is pulled towards the cities, one random city at a time,
until it settles into a closed curve that visits every city; reading the ring
in order gives the tour. It is the only solver that needs coordinates
(fit(X, coords=...)) and it never reads the cost matrix during its search — the
matrix prices the tour decoded after every epoch, which is what history_
records (best-so-far, SPEC R8).
Under a multi-trip budget the map ignores max_time_work (it warns, D6); the
returned tour is still split into trips and priced under the multi-trip
objective. Like every 2.0 solver it evaluates on a dense (n, n) cost matrix,
so the practical ceiling is about 20 000 nodes; coordinate-only fitting with
distances computed on the fly is a 2.1 item (D18).
>>> from skroute.metaheuristics import SOM
>>> from skroute.datasets import load_tsp
>>> wi = load_tsp("wi29")
>>> som = SOM(random_state=0).fit(wi.distance_matrix(), coords=wi.coords, labels=wi.labels)
>>> som.cost_ / wi.optimal_tour_length < 1.15
True
>>> som.stop_reason_ in {"converged", "max_iter"} and som.n_iter_ == len(som.history_)
True
SOM ¶
SOM(n_units: int | None = None, learning_rate: float = 0.8, lr_decay: float = 0.99997, radius: float | None = None, radius_decay: float = 0.9997, n_iter: int = 100000, random_state: int | Generator | None = None, verbose: int = 0)
Bases: BaseRouter
Self-organising map (Kohonen ring) for the TSP over planar coordinates.
A ring of n_units neurons is initialised at random inside the bounding box of the
normalised coordinates. Each sample presents one random city: the nearest neuron wins,
and every neuron moves towards the city by learning_rate times a Gaussian of its
wrapped ring distance to the winner (standard deviation radius); then
learning_rate *= lr_decay and radius *= radius_decay. Samples are grouped in
epochs of max(1, n_iter // 100) samples — the outer iterations — and after every
epoch the ring is decoded into a tour and priced with the problem's objective.
The solver needs coordinates (fit(X, coords=...), tag requires_coords): the
map is trained on the geometry only; the cost matrix X prices the decoded tours and
may be anything (asymmetric, road costs...).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
n_units
|
int or None
|
Neurons on the ring. |
None
|
learning_rate
|
float in (0, 1]
|
Initial step of the update towards the presented city; decays by |
0.8
|
lr_decay
|
float in (0, 1]
|
Multiplicative decay of |
0.99997
|
radius
|
float or None
|
Initial standard deviation (in ring positions) of the Gaussian neighbourhood around
the winner. |
None
|
radius_decay
|
float in (0, 1]
|
Multiplicative decay of |
0.9997
|
n_iter
|
int
|
Maximum number of samples (single-city presentations). They run in epochs of
|
100_000
|
random_state
|
(int, Generator or None)
|
Seed of the initial ring and of the city drawn at every sample. All randomness is
pre-drawn per epoch from |
None
|
verbose
|
int
|
|
0
|
Attributes:
| Name | Type | Description |
|---|---|---|
history_ |
ndarray of shape (n_iter_,), float64
|
Best-so-far cost after every epoch (monotone non-increasing; never the cost of the current ring, whose trace is not monotone — R8). |
n_iter_ |
int
|
Epochs actually run ( |
n_samples_ |
int
|
Samples actually drawn ( |
stop_reason_ |
{'converged', 'max_iter', 'callback'}
|
|
See `BaseRouter` for ``tour_``, ``route_``, ``trips_``, ``cost_`` |
|
|
and the other fitted attributes every solver shares. |
|
Notes
Algorithm (Angéniol et al., 1988; Kohonen, 1982). Coordinates are rescaled with
normalize_coords (the longer side of the bounding box spans
[0, 1], the aspect ratio is preserved, so radius and learning_rate mean the same
on every instance). Per sample: winner = argmin_j ||w_j - x||,
g_j = exp(-d(j, winner)^2 / (2 radius^2)) with d the wrapped ring distance, and
w_j += learning_rate * g_j * (x - w_j). Decoding maps every city to its closest neuron,
orders the cities by neuron index (ties by city index) and rotates the sequence to the depot.
The Gaussian's width is floored at 1e-3 ring positions (radius itself keeps decaying
and drives the stop rule): below that every neuron but the winner already gets g_j == 0
exactly in double precision, so the floor changes no finite result and only keeps extreme
radius/radius_decay values (radius_decay=0.5, radius=1e-200) from underflowing
radius**2 to 0.0 and turning the winner's update into 0/0.
Complexity. Every sample is O(n_units) vectorised numpy work; decoding an epoch is
O(n * n_units) in blocks of bounded memory, plus one O(n) evaluation. With the defaults
the radius reaches one ring position after roughly ln(0.8 n) / 0.0003 samples, so a
run usually stops by "converged" well before the 100 epochs.
Callback events (D30, D31). "start" has no tour (the ring has not been decoded yet),
with the extra keys radius, learning_rate and n_units; every epoch emits one
"iteration" whose tour is the epoch's decoded ring and best_tour the best epoch's,
with radius and learning_rate (after the epoch's decay), n_samples and ring —
an (n_units, 2) float64 copy of the neuron positions mapped back to the units of
problem.coords (the aspect-preserving min-max scaling inverted: w * span + min), so a
viewer can draw the ring closing on the cities. The ring is built only when a callback is set.
Supports: symmetric and asymmetric cost matrices (the matrix only prices the decoded
tours), coordinates required; stochastic; iterative. Multi-trip: not budget-aware (D6):
under max_time_work the map ignores the budget during its search and fit warns; the
result is still decoded into trips and priced under the multi-trip objective (history_
records that objective).
Ceiling. Like every 2.0 solver, SOM evaluates on a dense (n, n) cost matrix, so the
practical ceiling is about 20 000 nodes (3.2 GB of float64). The four bundled instances
above it (vm22775, sw24978, bm33708, ch71009) can be read and subsampled
(load_tsp(name, n_nodes=5000)) but not solved whole; coordinate-only fitting with
Euclidean costs computed on the fly (fit(None, coords=xy)) is a 2.1 item (D18).
MultiStart with threads gains little on this numpy-bound solver
(R9); prefer prefer="processes" there.
References
.. [1] T. Kohonen, "Self-organized formation of topologically correct feature maps", Biological Cybernetics 43 (1982) 59-69. .. [2] B. Angéniol, G. de la Croix Vaubois and J.-Y. Le Texier, "Self-organizing feature maps and the travelling salesman problem", Neural Networks 1 (1988) 289-293.
Examples:
>>> import numpy as np
>>> from skroute.metaheuristics import SOM
>>> from skroute.datasets import load_tsp
>>> wi = load_tsp("wi29") # Western Sahara, optimum 27603
>>> som = SOM(random_state=0).fit(wi.distance_matrix(), coords=wi.coords, labels=wi.labels)
>>> som.cost_ / wi.optimal_tour_length < 1.15 # the fast-tier tolerance of SPEC §6
True
>>> int(som.route_[0]) == int(som.route_[-1]) == int(som.depot_) == 1
True
>>> som.n_iter_ == len(som.history_) and som.stop_reason_ in {"converged", "max_iter"}
True
>>> bool(np.all(np.diff(som.history_) <= 0)) and som.n_samples_ <= som.n_iter
True
Without coordinates the solver refuses to fit:
>>> SOM().fit(wi.distance_matrix())
Traceback (most recent call last):
...
ValueError: SOM needs node coordinates: fit(X, coords=...)