5.2.3 Simulated Annealing

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124

/*

Approximately minimizes an energy function $E(s)$ over an arbitrary state space using simulated
annealing. This randomized heuristic is useful when the states may be discrete, the energy has many
local minima, and gradients or stronger structural properties are unavailable.

Starting from an initial state $s$, each iteration generates a random neighboring state $s'$. An
improving move is always accepted, while a worsening move is accepted with probability
$\exp(-(E(s') - E(s))/T)$ at temperature $T$. High temperatures encourage exploration across energy
barriers; as the temperature decreases, the search increasingly favors improvements. The best state
is tracked separately because the current state may later move to a worse one.

The geometric schedule multiplies the temperature by `cooling_rate` each round, performing
$k = \lceil \log(T_{end}/T_{start}) / \log(r) \rceil$ iterations. Temperature has the same scale as
energy, so scaling every energy requires scaling both temperature bounds by the same factor to
preserve acceptance probabilities.

- `anneal_min(initial, energy, rand_neighbor, rng, ...)` returns (`best_energy`, `best_state`) given
  an initial state `initial`, a callable `energy(state)` that returns the state's energy, and a
  callable `rand_neighbor(state, temp, rng)` that returns a randomly chosen nearby state. Optional
  parameters default to `temp_start = 1000`, `temp_end = 1e-6`, `cooling_rate = 0.995`, and
  `time_limit = 0`. `temp_end` must be positive and less than `temp_start`, while `cooling_rate`
  must be strictly between $0$ and $1$. A positive `time_limit`, in seconds, stops the search once
  that much wall-clock time has elapsed; the best state found so far is returned either way.

For continuous states, temperature can also control the neighbor step size. Slower cooling and
independent restarts often improve results, but simulated annealing never proves optimality.

This generic interface copies each proposed state and recomputes its energy. For large permutation
states, adapt the loop to mutate and revert moves in place and update the energy from the affected
terms instead.

Time Complexity:
- O(k(C_E + C_N)) per call, where $k$ is the number of temperature levels and $C_E$ and $C_N$ are
  the costs of evaluating the energy and generating a neighboring state. A positive `time_limit`
  caps $k$ at whatever the deadline allows.

Space Complexity:
- O(S) for the current, next, and best states, where $S$ is the size of one state.

*/

#include <cassert>
#include <chrono>
#include <cmath>
#include <random>
#include <utility>

template<typename State, typename Energy, typename Neighbor>
std::pair<double, State> anneal_min(
    State initial, Energy energy, Neighbor rand_neighbor, std::mt19937 &rng,
    double temp_start = 1000, double temp_end = 1e-6, double cooling_rate = 0.995,
    double time_limit = 0
) {
  assert(0 < temp_end && temp_end < temp_start);
  assert(0 < cooling_rate && cooling_rate < 1);
  std::uniform_real_distribution<double> unit(0.0, 1.0);
  State current = std::move(initial), best = current;
  double curr_energy = energy(current), best_energy = curr_energy;
  auto deadline = std::chrono::steady_clock::now() + std::chrono::duration<double>(time_limit);
  for (double temp = temp_start; temp > temp_end; temp *= cooling_rate) {
    if (time_limit > 0 && std::chrono::steady_clock::now() >= deadline) {
      break;  // The schedule is abandoned mid-cooling; the best state so far is still returned.
    }
    State next = rand_neighbor(current, temp, rng);
    double next_energy = energy(next);
    if (next_energy <= curr_energy || unit(rng) < std::exp((curr_energy - next_energy) / temp)) {
      current = std::move(next);
      curr_energy = next_energy;
      if (curr_energy < best_energy) {
        best = current;
        best_energy = curr_energy;
      }
    }
  }
  return {best_energy, std::move(best)};
}

/*** Example Usage ***/

#include <algorithm>
#include <vector>
using namespace std;

int main() {
  mt19937 rng(1234567);  // Fixed seed for reproducibility.

  // Escape the local minimum near x = 1 and find the lower minimum near x = -1.
  auto energy_1d = [](double x) { return (x * x - 1) * (x * x - 1) + 0.2 * x; };
  auto neighbor_1d = [](double x, double temp, mt19937 &gen) {
    uniform_real_distribution<double> move(-1.0, 1.0);
    return x + move(gen) * sqrt(temp);
  };
  auto [best_energy, best_x] = anneal_min(1.0, energy_1d, neighbor_1d, rng, 2.0, 1e-8);
  assert(best_energy < -0.2 && best_x < -1);

  // Approximate TSP with a permutation state and random 2-opt reversals.
  vector<pair<double, double>> points{{0, 0}, {1, 0}, {2, 0}, {2, 1},
                                      {2, 2}, {1, 2}, {0, 2}, {0, 1}};
  auto tour_length = [&](const vector<int> &tour) {
    double length = 0;
    for (int i = 0; i < static_cast<int>(tour.size()); i++) {
      auto [x1, y1] = points[tour[i]];
      auto [x2, y2] = points[tour[(i + 1) % tour.size()]];
      length += hypot(x1 - x2, y1 - y2);
    }
    return length;
  };
  auto reverse_segment = [](vector<int> tour, double, mt19937 &gen) {
    uniform_int_distribution<int> pick(0, static_cast<int>(tour.size()) - 1);
    int lo = pick(gen), hi = pick(gen);
    if (lo > hi) {
      swap(lo, hi);
    }
    reverse(tour.begin() + lo, tour.begin() + hi + 1);
    return tour;
  };
  vector<int> initial_tour{0, 4, 2, 6, 1, 5, 3, 7};
  auto [best_length, best_tour] =
      anneal_min(initial_tour, tour_length, reverse_segment, rng, 5.0, 1e-6, 0.995);
  assert(fabs(best_length - 8) < 1e-9);
  assert(fabs(tour_length(best_tour) - best_length) < 1e-9);
  return 0;
}