Skip to content

Simulated Annealing

General Principle

Simulated annealing is a probabilistic optimization technique inspired by the metallurgical process of annealing, where metals are heated and slowly cooled to reduce defects and reach a low-energy state.

The key idea is to allow the algorithm to:

  • Explore the solution space through random movements (like thermal fluctuations)
  • Exploit promising regions by moving toward better solutions
  • Gradually reduce randomness over time (cooling), converging to a good solution

This balance between exploration and exploitation helps escape local minima—a critical advantage over greedy algorithms like Lloyd’s method.

Mechanisms in K-means Context

In kmeanssa-ng, cluster centers are dynamic entities that move through the metric space using two complementary mechanisms:

Brownian Motion (Exploration)

Centers perform random walks in the metric space, characterized by:

  • Random direction selection
  • Distance traveled proportional to \(\sqrt{\Delta t}\) (diffusion scaling)
  • Allows escape from local minima through stochastic exploration

For a center \(c\) and time parameter \(\Delta t\):

\[ c \leftarrow c + \text{Brownian}(\Delta t) \]

Drift (Exploitation)

Centers are pulled toward the observations assigned to their cluster:

  • Each center drifts toward a randomly selected point in its cluster
  • Distance traveled: proportion \(\alpha\) of the geodesic distance
  • Reduces cluster energy by moving centers closer to their observations

For a center \(c\), target observation \(x\), and drift proportion \(\alpha \in [0,1]\):

\[ c \leftarrow c + \alpha \cdot (x - c) \]

(where the notation is geometric; on graphs this means moving along the geodesic path)

Temperature Schedule

The “temperature” controls the balance between exploration and exploitation over time. kmeanssa-ng advances an annealing clock with an inhomogeneous Poisson process of intensity \(\lambda(t) = \lambda_0 (1 + t)\) (the schedule of the companion paper). The \(n\)-th observation is processed at clock time

\[ T(n) = \sqrt{\,2\sum_{i=1}^{n} E_i + 1\,} - 1 \]

where \(E_i \sim \text{Exp}(\lambda_0)\) are the inter-arrival intervals.

Key properties:

  • Temperature decreases over time (cooling)
  • Controlled by parameter \(\lambda\) (intensity)
  • Stochastic schedule adds robustness

Running the Algorithm

The algorithm runs through a single method, run, which interleaves exploration and exploitation as it sweeps once through the (shuffled) observations. For each observation \(x_i\) it:

  1. Advances the annealing clock and performs Brownian motion on every center: \(c_j \leftarrow \text{move}(c_j, \Delta t)\)
  2. Finds the nearest center \(c^*\) to \(x_i\) and drifts it toward the observation

run takes an initialization strategy (how the centers start—e.g. KMeansPlusPlus), a robustification strategy (how the final centers are chosen from the tail of the sweep—e.g. MinimizeEnergy), and robust_prop, the fraction of the final observations used for that selection:

from kmeanssa_ng import generate_sbm, SimulatedAnnealing, KMeansPlusPlus, MinimizeEnergy
from kmeanssa_ng.quantum_graph.sampling import UniformNodeSampling

graph = generate_sbm(sizes=[25, 25], p=[[0.8, 0.1], [0.1, 0.8]], random_state=0)
points = graph.sample_points(150, strategy=UniformNodeSampling(random_state=0))

sa = SimulatedAnnealing(points, k=2, beta0=0.5, step_size=0.05, random_state=0)
centers = sa.run(KMeansPlusPlus(), MinimizeEnergy(), robust_prop=0.1)
print(f"Found {len(centers)} cluster centers")

Found 2 cluster centers

Passing an integer random_state (as above) makes a run fully reproducible.

Two robustification strategies are built in, trading quality for speed. MinimizeEnergy (the default) keeps the lowest-energy state in the collection window — the higher-quality choice, but it recomputes the energy repeatedly. On a quantum graph, MostFrequentNode instead returns the most frequently visited node, which is faster and places each center exactly on a node, at the cost of some selection accuracy. Prefer MinimizeEnergy unless the selection cost dominates your runtime.

Algorithm Parameters

The behavior of the simulated annealing algorithm is controlled by three main parameters. Understanding these parameters is crucial for achieving good clustering results.

lambda0: Observation-Clock Intensity

Controls the rate at which observations arrive, and therefore how much Brownian exploration happens between them.

Mathematical role: lambda0 is the scale of the inhomogeneous Poisson clock, of intensity \(\lambda(t) = \lambda_0 (1 + t)\); it sets the pace of the annealing schedule. It does not scale the Brownian steps — each micro-step has standard deviation \(\sqrt{\Delta t}\) (the step_size), independent of lambda0.

Practical effects:

  • Higher values (\(\lambda_0 \in [1.5, 3.0]\)):
  • Observations arrive faster, so the same number of observations spans a shorter annealing horizon
  • Less Brownian exploration between observations
  • Higher risk of getting trapped in a local minimum
  • Lower values (\(\lambda_0 \in [0.3, 0.8]\)):
  • A longer horizon, hence more Brownian exploration between observations
  • Better escape from local minima, at the cost of more computation
  • Default: \(\lambda_0 = 1.0\) provides a good balance for most use cases.

Choosing it for your problem: the annealing horizon grows like \(\sqrt{n_\text{obs} / \lambda_0}\), so with more observations you can afford a larger lambda0 (each observation already refines the centres). On larger graphs or with more clusters, keep lambda0 at or below \(1.0\): the longer horizon gives the centres the exploration they need to separate distant modes rather than collapsing onto the nearest one.

beta0: Drift Intensity

Controls how strongly centers are pulled toward observations.

Mathematical role: The drift proportion at time \(t\) is computed as: $\(\alpha(t) = \min(h \cdot \beta_0 \cdot \log(1 + t), 1)\)$ where \(h\) is the time interval. This controls the strength of attraction toward the nearest observation.

Practical effects:

  • Higher values (\(\beta_0 \in [2.0, 5.0]\)):
  • Stronger drift toward observations
  • Faster convergence
  • More exploitation of currently good positions
  • Risk of premature convergence
  • Lower values (\(\beta_0 \in [0.3, 0.8]\)):
  • Weaker drift
  • More time for exploration
  • Slower convergence
  • Better for difficult optimization landscapes
  • Default: \(\beta_0 \in [1.0, 2.0]\) works well for most problems

Interaction with lambda0: lambda0 sets how long the annealing runs (the horizon), while beta0 sets how hard each observation pulls its nearest centre. A strong drift (beta0) on a short horizon (large lambda0) converges fast but can lock the centres in before they have explored — the classic premature-convergence failure. If clusters come out merged or a centre is stranded, lower beta0 or lambda0 (lengthening the horizon) rather than raising them.

step_size: Time Discretization

Controls the temporal resolution of the stochastic process simulation.

Mathematical role: Euler discretization step \(\Delta t\) for solving the stochastic differential equation (SDE). Smaller values provide more accurate simulation of the continuous process.

Practical effects:

  • Smaller values (\(\Delta t \in [0.001, 0.01]\)):
  • More accurate simulation of the theoretical process
  • Slower computation (more steps needed)
  • Better numerical stability
  • Larger values (\(\Delta t \in [0.05, 0.1]\)):
  • Faster computation
  • Less accurate approximation
  • Risk of numerical instability
  • Default: \(\Delta t = 0.01\) provides a good accuracy/speed tradeoff

Rule of thumb: Choose step_size much smaller than the characteristic time scale of your Poisson process (approximately \(1/\lambda_0\)).

Checking it is small enough: a Brownian micro-step moves a centre by about \(\sqrt{\Delta t}\). If that is comparable to the spacing between clusters, a centre can jump across a cluster boundary in a single step and the discretisation is too coarse. Halve step_size until the clustering stops changing; as a guide, keep \(\sqrt{\Delta t}\) well below the smallest inter-cluster distance in your space.

Choosing Parameter Combinations

Different parameter combinations suit different optimization scenarios:

Quick Convergence (when you trust your initialization)

sa = SimulatedAnnealing(
    points, k=5,
    lambda0=0.5,    # Reduced exploration
    beta0=3.0,      # Strong drift
    step_size=0.01
)

Thorough Search (complex problems, many local minima)

sa = SimulatedAnnealing(
    points, k=5,
    lambda0=2.0,    # More exploration
    beta0=1.0,      # Gentler drift
    step_size=0.01
)

Balanced (default, works well in most cases)

sa = SimulatedAnnealing(
    points, k=5,
    lambda0=1.0,    # Moderate exploration
    beta0=2.0,      # Moderate drift
    step_size=0.01
)

Theoretical Foundation

For the derivation of the annealing dynamics and its convergence analysis, see the companion paper: C. Brécheteau, I. Gavra, N. Klutchnikoff, Online k-means Clustering on Metric Graphs and Geodesic Spaces (in preparation).

Monitoring Convergence

Pass record_energy=True to run to record the \(k\)-means energy after every observation (plus the initial state). The trajectory is exposed as energy_history, with the matching annealing times as time_history. It is the raw exploration path of the centers—it rises and falls as the anneal explores—and the robustification returns its lowest-energy configuration. Plotting it is the usual way to see how the anneal progresses:

sa = SimulatedAnnealing(points, k=2, beta0=0.5, step_size=0.05, random_state=0)
sa.run(KMeansPlusPlus(), MinimizeEnergy(), robust_prop=0.1, record_energy=True)
history = sa.energy_history          # energy after each observation
print(f"recorded {len(history)} energies; lowest reached = {history.min():.3f}")

recorded 151 energies; lowest reached = 1.224

Recording does not affect the result and is off by default, so the energy is only recomputed when you ask for it.

Robustification

To improve stability, kmeanssa-ng uses robustification: instead of returning the final centers, it averages results from the last \(p\%\) of iterations (default: 10%). This reduces sensitivity to late-iteration fluctuations.