4.7.4 Shortest Hamiltonian Path and Cycle (Held-Karp)

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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

/*

Given a complete, weighted, directed graph, find a minimum-length Hamiltonian path or cycle. A
Hamiltonian path visits each node exactly once. A Hamiltonian cycle additionally returns to its
starting node; this cycle version is the traveling salesman problem (TSP).

Both functions use the Held-Karp subset dynamic program. Let $dp(S, i)$ be the shortest path that
visits exactly the nodes in $S$ and ends at $i$. Its final edge must come from some
$j \in S \setminus \{i\}$, giving the recurrence
$dp(S, i) = \min_j(dp(S \setminus \{i\}, j) + w(j, i))$. Paths initialize every singleton set;
cycles initialize only node $0$ and consider only subsets containing it. After processing the full
set, the algorithm chooses the best final node, including its return edge to node $0$ for a cycle.
Reconstruction repeats the same minimization backward instead of storing a separate parent table.

- `shortest_hamiltonian_path()` populates `path` and returns the minimum Hamiltonian path length for
  a global, pre-populated adjacency matrix `adj`.
- `shortest_hamiltonian_cycle()` fixes the start at node 0, populates `path`, and returns the
  minimum Hamiltonian cycle length for a global, pre-populated adjacency matrix `adj`. The path
  lists every node once; the closing edge back to node 0 is implicit.

Since this implementation uses bitmasks with signed 32-bit integers, the maximum number of nodes
must be less than $31$.

Time Complexity:
- O(2^n * n^2) per call, where $n$ is the number of nodes.

Space Complexity:
- O(n^2) for storage of the graph, where $n$ is the number of nodes.
- O(2^n * n) auxiliary.

*/

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <vector>

const int64_t INF = INT64_MAX / 4;
std::vector<std::vector<int64_t>> adj;
std::vector<int> path;

int64_t shortest_hamiltonian_path() {
  int n = static_cast<int>(adj.size());
  assert(1 <= n && n < 31);
  int max_mask = (1 << n) - 1;
  std::vector<std::vector<int64_t>> dp(max_mask + 1, std::vector<int64_t>(n, INF));
  path.assign(n, 0);
  for (int i = 0; i < n; i++) {
    dp[1 << i][i] = 0;
  }
  for (int mask = 1; mask <= max_mask; mask++) {
    for (int i = 0; i < n; i++) {
      if ((mask & (1 << i)) != 0) {
        for (int j = 0; j < n; j++) {
          if ((mask & (1 << j)) != 0 && dp[mask ^ (1 << i)][j] != INF) {
            int64_t candidate = dp[mask ^ (1 << i)][j] + adj[j][i];  // Overflow warning.
            dp[mask][i] = std::min(dp[mask][i], candidate);
          }
        }
      }
    }
  }
  int64_t res = INF;
  for (int i = 0; i < n; i++) {
    res = std::min(res, dp[max_mask][i]);
  }
  // Optional: reconstruct one shortest Hamiltonian path.
  int mask = max_mask, old = -1;
  for (int i = n - 1; i >= 0; i--) {
    int best = -1;
    for (int j = 0; j < n; j++) {
      if ((mask & (1 << j)) != 0 &&
          (best == -1 || dp[mask][best] + (old == -1 ? 0 : adj[best][old]) >
                             dp[mask][j] + (old == -1 ? 0 : adj[j][old]))) {
        best = j;
      }
    }
    path[i] = best;
    mask ^= 1 << best;
    old = best;
  }
  return res;
}

int64_t shortest_hamiltonian_cycle() {
  int n = static_cast<int>(adj.size());
  assert(1 <= n && n < 31);
  if (n == 1) {
    path = {0};
    return 0;
  }
  int max_mask = (1 << n) - 1;
  std::vector<std::vector<int64_t>> dp(max_mask + 1, std::vector<int64_t>(n, INF));
  path.assign(n, 0);
  dp[1][0] = 0;
  for (int mask = 1; mask <= max_mask; mask += 2) {
    for (int i = 1; i < n; i++) {
      if ((mask & (1 << i)) != 0) {
        for (int j = 0; j < n; j++) {
          if ((mask & (1 << j)) != 0 && dp[mask ^ (1 << i)][j] != INF) {
            int64_t candidate = dp[mask ^ (1 << i)][j] + adj[j][i];  // Overflow warning.
            dp[mask][i] = std::min(dp[mask][i], candidate);
          }
        }
      }
    }
  }
  int64_t res = INF;
  for (int i = 1; i < n; i++) {
    res = std::min(res, dp[max_mask][i] + adj[i][0]);  // Overflow warning.
  }
  // Optional: reconstruct one shortest Hamiltonian cycle.
  int mask = max_mask, old = 0;
  for (int i = n - 1; i >= 1; i--) {
    int best = -1;
    for (int j = 1; j < n; j++) {
      if ((mask & (1 << j)) != 0 &&
          (best == -1 || dp[mask][best] + adj[best][old] > dp[mask][j] + adj[j][old])) {
        best = j;
      }
    }
    path[i] = best;
    mask ^= 1 << best;
    old = best;
  }
  return res;
}

/*** Example Usage ***/

#include <cassert>
using namespace std;

int main() {
  // Complete directed graph where 0->1->2->0 costs 1 + 2 + 3.
  // The reverse cycle is intentionally more expensive.
  //    +--> 0 <-------+
  //    |    | \       |
  // w=7| w=1|   \w=1  |w=3
  //    |    |     \   |
  //    |    v  w=2  v |
  //    +--- 1 ------> 2
  //          ^_______/
  //             w=5
  int nodes = 3;
  adj.assign(nodes, vector<int64_t>(nodes));
  adj[0][1] = 1;
  adj[0][2] = 1;
  adj[1][0] = 7;
  adj[1][2] = 2;
  adj[2][0] = 3;
  adj[2][1] = 5;
  assert(shortest_hamiltonian_path() == 3);
  assert((path == vector<int>{0, 1, 2}));
  assert(shortest_hamiltonian_cycle() == 6);
  assert((path == vector<int>{0, 1, 2}));
  return 0;
}