1 // Copyright 2010-2021 Google LLC
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 // http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13
14 // An implementation of the Held-Karp symmetric Traveling Salesman (TSP) lower
15 // bound algorithm, inspired by "Estimating the Held-Karp lower bound for the
16 // geometric TSP" by Christine L. Valenzuela and Antonia J. Jones, European
17 // Journal of Operational Research, Volume 102, Issue 1, 1 October 1997,
18 // Pages 157-175.
19 //
20 // The idea is to compute minimum 1-trees to evaluate a lower bound to the
21 // corresponding TSP. A minimum 1-tree is a minimum spanning tree on all nodes
22 // but one, to which are added the two shortest edges from the left-out node to
23 // the nodes of the spanning tree. The sum of the cost of the edges of the
24 // minimum 1-tree is a lower bound to the cost of the TSP.
25 // In order to improve (increase) this lower bound, the idea is to add weights
26 // to each nodes, weights which are added to the cost function used when
27 // computing the 1-tree. If weight[i] is the weight of node i, the cost function
28 // therefore becomes weighed_cost(i,j) = cost(i,j) + weight[i] + weight[j]. One
29 // can see that w = weighed_cost(minimum 1-tree) - Sum(2 * weight[i])
30 // = cost(minimum 1-tree) + Sum(weight[i] * (degree[i] - 2))
31 // is a valid lower bound to the TSP:
32 // 1) let T be the set of 1-trees on the nodes;
33 // 2) let U be the set of tours on the nodes; U is a subset of T (tours are
34 // 1-trees with all degrees equal to 2), therefore:
35 // min(t in T) Cost(t) <= min(t in U) Cost(t)
36 // and
37 // min(t in T) WeighedCost(t) <= min(t in U) WeighedCost(t)
38 // 3) weighed_cost(i,j) = cost(i,j) + weight[i] + weight[j], therefore:
39 // for all t in T, WeighedCost(t) = Cost(t) + Sum(weight[i] * degree[i])
40 // and
41 // for all i in U, WeighedCost(t) = Cost(t) + Sum(weight[i] * 2)
42 // 4) let t* in U s.t. WeighedCost(t*) = min(t in U) WeighedCost(t), therefore:
43 // min(t in T) (Cost(t) + Sum(weight[i] * degree[i]))
44 // <= Cost(t*) + Sum(weight[i] * 2)
45 // and
46 // min(t in T) (Cost(t) + Sum(weight[i] * (degree[i] - 2))) <= Cost(t*)
47 // and
48 // cost(minimum 1-tree) + Sum(weight[i] * (degree[i] - 2)) <= Cost(t*)
49 // and
50 // w <= Cost(t*)
51 // 5) because t* is also the tour minimizing Cost(t) with t in U (weights do not
52 // affect the optimality of a tour), Cost(t*) is the cost of the optimal
53 // solution to the TSP and w is a lower bound to this cost.
54 //
55 // The best lower bound is the one for which weights maximize w. Intuitively as
56 // degrees get closer to 2 the minimum 1-trees gets closer to a tour.
57 //
58 // At each iteration m, weights are therefore updated as follows:
59 // weight(m+1)[i] = weight(m)[i] + step(m) * (degree(m)[i] - 2)
60 // where degree(m)[i] is the degree of node i in the 1-tree at iteration i,
61 // step(m) is a subgradient optimization step.
62 //
63 // This implementation uses two variants of Held-Karp's initial subgradient
64 // optimization iterative estimation approach described in "The
65 // traveling-salesman problem and minimum spanning trees: Part I and II", by
66 // Michael Held and Richard M. Karp, Operations Research Vol. 18,
67 // No. 6 (Nov. - Dec., 1970), pp. 1138-1162 and Mathematical Programming (1971).
68 //
69 // The first variant comes from Volgenant, T., and Jonker, R. (1982), "A branch
70 // and bound algorithm for the symmetric traveling salesman problem based on the
71 // 1-tree relaxation", European Journal of Operational Research. 9:83-89.".
72 // It suggests using
73 // step(m) = (1.0 * (m - 1) * (2 * M - 5) / (2 * (M - 1))) * step1
74 // - (m - 2) * step1
75 // + (0.5 * (m - 1) * (m - 2) / ((M - 1) * (M - 2))) * step1
76 // where M is the maximum number of iterations and step1 is initially set to
77 // L / (2 * number of nodes), where L is the un-weighed cost of the 1-tree;
78 // step1 is updated each time a better w is found. The intuition is to have a
79 // positive decreasing step which is equal to 0 after M iterations; Volgenant
80 // and Jonker suggest that:
81 // step(m) - 2 * step(m-1) + t(m-2) = constant,
82 // step(M) = 0
83 // and
84 // step(1) - step(2) = 3 * (step(M-1) - step(M)).
85 // The step(m) formula above derives from this recursive formulation.
86 // This is the default algorithm used in this implementation.
87 //
88 // The second variant comes from Held, M., Wolfe, P., and Crowder, H. P. (1974),
89 // "Validation of subgradient optimization", Mathematical Programming 6:62-88.
90 // It derives from the original Held-Karp formulation:
91 // step(m) = lambda(m) * (wlb - w(m)) / Sum((degree[i] - 2)^2),
92 // where wlb is a lower bound to max(w(m)) and lambda(m) in [0, 2].
93 // Help-Karp prove that
94 // if w(m') > w(m) and 0 < step < 2 * (w(m') - w(m))/norm(degree(m) - 2)^2,
95 // then weight(m+1) is closer to w' than w from which they derive the above
96 // formula.
97 // Held-Wolfe-Crowder show that using an overestimate UB is as effective as
98 // using the underestimate wlb while UB is easier to compute. The resulting
99 // formula is:
100 // step(m) = lambda(m) * (UB - w(m)) / Sum((degree[i] - 2)^2),
101 // where UB is an upper bound to the TSP (here computed with the Christofides
102 // algorithm), and lambda(m) in [0, 2] initially set to 2. Held-Wolfe-Crowder
103 // suggest running the algorithm for M = 2 * number of nodes iterations, then
104 // dividing lambda and M by 2 until M is small enough (less than 2 in this
105 // implementation).
106 //
107 // To speed up the computation, minimum spanning trees are actually computed on
108 // a graph limited to the nearest neighbors of each node. Valenzuela-Jones 1997
109 // experiments have shown that this does not harm the lower bound computation
110 // significantly. At the end of the algorithm a last iteration is run on the
111 // complete graph to ensure the bound is correct (the cost of a minimum 1-tree
112 // on a partial graph is an upper bound to the one on a complete graph).
113 //
114 // Usage:
115 // std::function<int64_t(int,int)> cost_function =...;
116 // const double lower_bound =
117 // ComputeOneTreeLowerBound(number_of_nodes, cost_function);
118 // where number_of_nodes is the number of nodes in the TSP and cost_function
119 // is a function returning the cost between two nodes.
120
121 #ifndef OR_TOOLS_GRAPH_ONE_TREE_LOWER_BOUND_H_
122 #define OR_TOOLS_GRAPH_ONE_TREE_LOWER_BOUND_H_
123
124 #include <math.h>
125
126 #include <cstdint>
127 #include <limits>
128 #include <set>
129
130 #include "ortools/base/integral_types.h"
131 #include "ortools/graph/christofides.h"
132 #include "ortools/graph/minimum_spanning_tree.h"
133
134 namespace operations_research {
135
136 // Implementation of algorithms computing Held-Karp bounds. They have to provide
137 // the following methods:
138 // - bool Next(): returns false when the algorithm must stop;
139 // - double GetStep(): returns the current step computed by the algorithm;
140 // - void OnOneTree(CostType one_tree_cost,
141 // double w,
142 // const std::vector<int>& degrees):
143 // called each time a new minimum 1-tree is computed;
144 // - one_tree_cost: the un-weighed cost of the 1-tree,
145 // - w the current value of w,
146 // - degrees: the degree of nodes in the 1-tree.
147 // - OnNewWMax(CostType one_tree_cost): called when a better value of w is
148 // found, one_tree_cost being the un-weighed cost of the corresponding
149 // minimum 1-tree.
150
151 // Implementation of the Volgenant Jonker algorithm (see the comments at the
152 // head of the file for explanations).
153 template <typename CostType>
154 class VolgenantJonkerEvaluator {
155 public:
VolgenantJonkerEvaluator(int number_of_nodes,int max_iterations)156 VolgenantJonkerEvaluator(int number_of_nodes, int max_iterations)
157 : step1_initialized_(false),
158 step1_(0),
159 iteration_(0),
160 max_iterations_(max_iterations > 0 ? max_iterations
161 : MaxIterations(number_of_nodes)),
162 number_of_nodes_(number_of_nodes) {}
163
Next()164 bool Next() { return iteration_++ < max_iterations_; }
165
GetStep()166 double GetStep() const {
167 return (1.0 * (iteration_ - 1) * (2 * max_iterations_ - 5) /
168 (2 * (max_iterations_ - 1))) *
169 step1_ -
170 (iteration_ - 2) * step1_ +
171 (0.5 * (iteration_ - 1) * (iteration_ - 2) /
172 ((max_iterations_ - 1) * (max_iterations_ - 2))) *
173 step1_;
174 }
175
OnOneTree(CostType one_tree_cost,double w,const std::vector<int> & degrees)176 void OnOneTree(CostType one_tree_cost, double w,
177 const std::vector<int>& degrees) {
178 if (!step1_initialized_) {
179 step1_initialized_ = true;
180 UpdateStep(one_tree_cost);
181 }
182 }
183
OnNewWMax(CostType one_tree_cost)184 void OnNewWMax(CostType one_tree_cost) { UpdateStep(one_tree_cost); }
185
186 private:
187 // Automatic computation of the number of iterations based on empirical
188 // results given in Valenzuela-Jones 1997.
MaxIterations(int number_of_nodes)189 static int MaxIterations(int number_of_nodes) {
190 return static_cast<int>(28 * std::pow(number_of_nodes, 0.62));
191 }
192
UpdateStep(CostType one_tree_cost)193 void UpdateStep(CostType one_tree_cost) {
194 step1_ = one_tree_cost / (2 * number_of_nodes_);
195 }
196
197 bool step1_initialized_;
198 double step1_;
199 int iteration_;
200 const int max_iterations_;
201 const int number_of_nodes_;
202 };
203
204 // Implementation of the Held-Wolfe-Crowder algorithm (see the comments at the
205 // head of the file for explanations).
206 template <typename CostType, typename CostFunction>
207 class HeldWolfeCrowderEvaluator {
208 public:
HeldWolfeCrowderEvaluator(int number_of_nodes,const CostFunction & cost)209 HeldWolfeCrowderEvaluator(int number_of_nodes, const CostFunction& cost)
210 : iteration_(0),
211 number_of_iterations_(2 * number_of_nodes),
212 upper_bound_(0),
213 lambda_(2.0),
214 step_(0) {
215 // TODO(user): Improve upper bound with some local search; tighter upper
216 // bounds lead to faster convergence.
217 ChristofidesPathSolver<CostType, int64_t, int, CostFunction> solver(
218 number_of_nodes, cost);
219 upper_bound_ = solver.TravelingSalesmanCost();
220 }
221
Next()222 bool Next() {
223 const int min_iterations = 2;
224 if (iteration_ >= number_of_iterations_) {
225 number_of_iterations_ /= 2;
226 if (number_of_iterations_ < min_iterations) return false;
227 iteration_ = 0;
228 lambda_ /= 2;
229 } else {
230 ++iteration_;
231 }
232 return true;
233 }
234
GetStep()235 double GetStep() const { return step_; }
236
OnOneTree(CostType one_tree_cost,double w,const std::vector<int> & degrees)237 void OnOneTree(CostType one_tree_cost, double w,
238 const std::vector<int>& degrees) {
239 double norm = 0;
240 for (int degree : degrees) {
241 const double delta = degree - 2;
242 norm += delta * delta;
243 }
244 step_ = lambda_ * (upper_bound_ - w) / norm;
245 }
246
OnNewWMax(CostType one_tree_cost)247 void OnNewWMax(CostType one_tree_cost) {}
248
249 private:
250 int iteration_;
251 int number_of_iterations_;
252 CostType upper_bound_;
253 double lambda_;
254 double step_;
255 };
256
257 // Computes the nearest neighbors of each node for the given cost function.
258 // The ith element of the returned vector contains the indices of the nearest
259 // nodes to node i. Note that these indices contain the number_of_neighbors
260 // nearest neighbors as well as all the nodes for which i is a nearest
261 // neighbor.
262 template <typename CostFunction>
NearestNeighbors(int number_of_nodes,int number_of_neighbors,const CostFunction & cost)263 std::set<std::pair<int, int>> NearestNeighbors(int number_of_nodes,
264 int number_of_neighbors,
265 const CostFunction& cost) {
266 using CostType = decltype(cost(0, 0));
267 std::set<std::pair<int, int>> nearest;
268 for (int i = 0; i < number_of_nodes; ++i) {
269 std::vector<std::pair<CostType, int>> neighbors;
270 neighbors.reserve(number_of_nodes - 1);
271 for (int j = 0; j < number_of_nodes; ++j) {
272 if (i != j) {
273 neighbors.emplace_back(cost(i, j), j);
274 }
275 }
276 int size = neighbors.size();
277 if (number_of_neighbors < size) {
278 std::nth_element(neighbors.begin(),
279 neighbors.begin() + number_of_neighbors - 1,
280 neighbors.end());
281 size = number_of_neighbors;
282 }
283 for (int j = 0; j < size; ++j) {
284 nearest.insert({i, neighbors[j].second});
285 nearest.insert({neighbors[j].second, i});
286 }
287 }
288 return nearest;
289 }
290
291 // Let G be the complete graph on nodes in [0, number_of_nodes - 1]. Adds arcs
292 // from the minimum spanning tree of G to the arcs set argument.
293 template <typename CostFunction>
AddArcsFromMinimumSpanningTree(int number_of_nodes,const CostFunction & cost,std::set<std::pair<int,int>> * arcs)294 void AddArcsFromMinimumSpanningTree(int number_of_nodes,
295 const CostFunction& cost,
296 std::set<std::pair<int, int>>* arcs) {
297 util::CompleteGraph<int, int> graph(number_of_nodes);
298 const std::vector<int> mst =
299 BuildPrimMinimumSpanningTree(graph, [&cost, &graph](int arc) {
300 return cost(graph.Tail(arc), graph.Head(arc));
301 });
302 for (int arc : mst) {
303 arcs->insert({graph.Tail(arc), graph.Head(arc)});
304 arcs->insert({graph.Head(arc), graph.Tail(arc)});
305 }
306 }
307
308 // Returns the index of the node in graph which minimizes cost(node, source)
309 // with the constraint that accept(node) is true.
310 template <typename CostFunction, typename GraphType, typename AcceptFunction>
GetNodeMinimizingEdgeCostToSource(const GraphType & graph,int source,const CostFunction & cost,AcceptFunction accept)311 int GetNodeMinimizingEdgeCostToSource(const GraphType& graph, int source,
312 const CostFunction& cost,
313 AcceptFunction accept) {
314 int best_node = -1;
315 double best_edge_cost = 0;
316 for (const auto node : graph.AllNodes()) {
317 if (accept(node)) {
318 const double edge_cost = cost(node, source);
319 if (best_node == -1 || edge_cost < best_edge_cost) {
320 best_node = node;
321 best_edge_cost = edge_cost;
322 }
323 }
324 }
325 return best_node;
326 }
327
328 // Computes a 1-tree for the given graph, cost function and node weights.
329 // Returns the degree of each node in the 1-tree and the un-weighed cost of the
330 // 1-tree.
331 template <typename CostFunction, typename GraphType, typename CostType>
ComputeOneTree(const GraphType & graph,const CostFunction & cost,const std::vector<double> & weights,const std::vector<int> & sorted_arcs,CostType * one_tree_cost)332 std::vector<int> ComputeOneTree(const GraphType& graph,
333 const CostFunction& cost,
334 const std::vector<double>& weights,
335 const std::vector<int>& sorted_arcs,
336 CostType* one_tree_cost) {
337 const auto weighed_cost = [&cost, &weights](int from, int to) {
338 return cost(from, to) + weights[from] + weights[to];
339 };
340 // Compute MST on graph.
341 std::vector<int> mst;
342 if (!sorted_arcs.empty()) {
343 mst = BuildKruskalMinimumSpanningTreeFromSortedArcs<GraphType>(graph,
344 sorted_arcs);
345 } else {
346 mst = BuildPrimMinimumSpanningTree<GraphType>(
347 graph, [&weighed_cost, &graph](int arc) {
348 return weighed_cost(graph.Tail(arc), graph.Head(arc));
349 });
350 }
351 std::vector<int> degrees(graph.num_nodes() + 1, 0);
352 *one_tree_cost = 0;
353 for (int arc : mst) {
354 degrees[graph.Head(arc)]++;
355 degrees[graph.Tail(arc)]++;
356 *one_tree_cost += cost(graph.Tail(arc), graph.Head(arc));
357 }
358 // Add 2 cheapest edges from the nodes in the graph to the extra node not in
359 // the graph.
360 const int extra_node = graph.num_nodes();
361 const auto update_one_tree = [extra_node, one_tree_cost, °rees,
362 &cost](int node) {
363 *one_tree_cost += cost(node, extra_node);
364 degrees.back()++;
365 degrees[node]++;
366 };
367 const int node = GetNodeMinimizingEdgeCostToSource(
368 graph, extra_node, weighed_cost,
369 [extra_node](int n) { return n != extra_node; });
370 update_one_tree(node);
371 update_one_tree(GetNodeMinimizingEdgeCostToSource(
372 graph, extra_node, weighed_cost,
373 [extra_node, node](int n) { return n != extra_node && n != node; }));
374 return degrees;
375 }
376
377 // Computes the lower bound of a TSP using a given subgradient algorithm.
378 template <typename CostFunction, typename Algorithm>
ComputeOneTreeLowerBoundWithAlgorithm(int number_of_nodes,int nearest_neighbors,const CostFunction & cost,Algorithm * algorithm)379 double ComputeOneTreeLowerBoundWithAlgorithm(int number_of_nodes,
380 int nearest_neighbors,
381 const CostFunction& cost,
382 Algorithm* algorithm) {
383 if (number_of_nodes < 2) return 0;
384 if (number_of_nodes == 2) return cost(0, 1) + cost(1, 0);
385 using CostType = decltype(cost(0, 0));
386 auto nearest = NearestNeighbors(number_of_nodes - 1, nearest_neighbors, cost);
387 // Ensure nearest arcs result in a connected graph by adding arcs from the
388 // minimum spanning tree; this will add arcs which are likely to be "good"
389 // 1-tree arcs.
390 AddArcsFromMinimumSpanningTree(number_of_nodes - 1, cost, &nearest);
391 util::ListGraph<int, int> graph(number_of_nodes - 1, nearest.size());
392 for (const auto& arc : nearest) {
393 graph.AddArc(arc.first, arc.second);
394 }
395 std::vector<double> weights(number_of_nodes, 0);
396 std::vector<double> best_weights(number_of_nodes, 0);
397 double max_w = -std::numeric_limits<double>::infinity();
398 double w = 0;
399 // Iteratively compute lower bound using a partial graph.
400 while (algorithm->Next()) {
401 CostType one_tree_cost = 0;
402 const std::vector<int> degrees =
403 ComputeOneTree(graph, cost, weights, {}, &one_tree_cost);
404 algorithm->OnOneTree(one_tree_cost, w, degrees);
405 w = one_tree_cost;
406 for (int j = 0; j < number_of_nodes; ++j) {
407 w += weights[j] * (degrees[j] - 2);
408 }
409 if (w > max_w) {
410 max_w = w;
411 best_weights = weights;
412 algorithm->OnNewWMax(one_tree_cost);
413 }
414 const double step = algorithm->GetStep();
415 for (int j = 0; j < number_of_nodes; ++j) {
416 weights[j] += step * (degrees[j] - 2);
417 }
418 }
419 // Compute lower bound using the complete graph on the best weights. This is
420 // necessary as the MSTs computed on nearest neighbors is not guaranteed to
421 // lead to a lower bound.
422 util::CompleteGraph<int, int> complete_graph(number_of_nodes - 1);
423 CostType one_tree_cost = 0;
424 // TODO(user): We are not caching here since this would take O(n^2) memory;
425 // however the Kruskal algorithm will expand all arcs also consuming O(n^2)
426 // memory; investigate alternatives to expanding all arcs (Prim's algorithm).
427 const std::vector<int> degrees =
428 ComputeOneTree(complete_graph, cost, best_weights, {}, &one_tree_cost);
429 w = one_tree_cost;
430 for (int j = 0; j < number_of_nodes; ++j) {
431 w += best_weights[j] * (degrees[j] - 2);
432 }
433 return w;
434 }
435
436 // Parameters to configure the computation of the TSP lower bound.
437 struct TravelingSalesmanLowerBoundParameters {
438 enum Algorithm {
439 VolgenantJonker,
440 HeldWolfeCrowder,
441 };
442 // Subgradient algorithm to use to compute the TSP lower bound.
443 Algorithm algorithm = VolgenantJonker;
444 // Number of iterations to use in the Volgenant-Jonker algorithm. Overrides
445 // automatic iteration computation if positive.
446 int volgenant_jonker_iterations = 0;
447 // Number of nearest neighbors to consider in the miminum spanning trees.
448 int nearest_neighbors = 40;
449 };
450
451 // Computes the lower bound of a TSP using given parameters.
452 template <typename CostFunction>
ComputeOneTreeLowerBoundWithParameters(int number_of_nodes,const CostFunction & cost,const TravelingSalesmanLowerBoundParameters & parameters)453 double ComputeOneTreeLowerBoundWithParameters(
454 int number_of_nodes, const CostFunction& cost,
455 const TravelingSalesmanLowerBoundParameters& parameters) {
456 using CostType = decltype(cost(0, 0));
457 switch (parameters.algorithm) {
458 case TravelingSalesmanLowerBoundParameters::VolgenantJonker: {
459 VolgenantJonkerEvaluator<CostType> algorithm(
460 number_of_nodes, parameters.volgenant_jonker_iterations);
461 return ComputeOneTreeLowerBoundWithAlgorithm(
462 number_of_nodes, parameters.nearest_neighbors, cost, &algorithm);
463 break;
464 }
465 case TravelingSalesmanLowerBoundParameters::HeldWolfeCrowder: {
466 HeldWolfeCrowderEvaluator<CostType, CostFunction> algorithm(
467 number_of_nodes, cost);
468 return ComputeOneTreeLowerBoundWithAlgorithm(
469 number_of_nodes, parameters.nearest_neighbors, cost, &algorithm);
470 }
471 default:
472 LOG(ERROR) << "Unsupported algorithm: " << parameters.algorithm;
473 return 0;
474 }
475 }
476
477 // Computes the lower bound of a TSP using default parameters (Volgenant-Jonker
478 // algorithm, 200 iterations and 40 nearest neighbors) which have turned out to
479 // give good results on the TSPLIB.
480 template <typename CostFunction>
ComputeOneTreeLowerBound(int number_of_nodes,const CostFunction & cost)481 double ComputeOneTreeLowerBound(int number_of_nodes, const CostFunction& cost) {
482 TravelingSalesmanLowerBoundParameters parameters;
483 return ComputeOneTreeLowerBoundWithParameters(number_of_nodes, cost,
484 parameters);
485 }
486
487 } // namespace operations_research
488
489 #endif // OR_TOOLS_GRAPH_ONE_TREE_LOWER_BOUND_H_
490