淺談旅行推銷員問題 (Travelling Salesman Problem) 搜索加速

contents

  1. 1. Task Description
  2. 2. Subtask
  3. 3. Input Format
  4. 4. Output Format
  5. 5. Sample Input 1
  6. 6. Sample Output 1
  7. 7. Solution
    1. 7.1. 初學第一步-窮舉
    2. 7.2. 踏出第一步-動態規劃
    3. 7.3. 更進一步
      1. 7.3.1. 計算基礎花費
      2. 7.3.2. 優化一夜
      3. 7.3.3. 優化二夜
      4. 7.3.4. 優化三夜
      5. 7.3.5. 優化四夜

學弟們修課〈演算法設計與分析〉內容中的一環,看他們的講義進行加速比較,也許有人會認為我太閒,但這只是追求真理。看看維基百科上,就有那上萬點的極佳近似解算法,而接下來要說的相當於玩具等級。

Task Description

給一組城市和每一個城市之間的單向距離,請找出從起點城市出發,經過所有城市且不重複回到起點城市的最短距離和。

例如,走訪順序為 $1-2-4-3-1$. 距離總和為 $10+25+30+15 = 80$.

Subtask

$N \le 100$

Input Format

只有一組測資,每組第一行有一個整數 $N$ ,表示有 $N$ 座城市,接下來會有 $N$ 行,每行上有 $N$ 個非負整數 $D_{i, j}$ 表示節點 $i$ 到節點 $j$ 的距離,如果 $D_{i, j} = 0$ ,視同沒有邊。

你可以假設每一條邊可能是單向的,保證至少一條 TSP 路徑。

  • $3 < N \leq 100$
  • $0 < \textit{distance between any pair of cities} \le 10000$

Output Format

輸出一行最小距離和。

Sample Input 1

1
2
3
4
5
4
0 10 15 20
10 0 35 25
15 35 0 30
20 25 30 0

Sample Output 1

1
80

Solution

初學第一步-窮舉

從最直觀的想法中,我們可以使用窮舉 $O(N!)$ 解決,再搭配簡單的剪枝 (加上當前點 $u$ 返回起點的最短路徑長是否大於已知解)。當 $N \ge 12$ 時,理論運算次數達 4 億之多,運算時間相當於 1 秒左右 (粗估一般 CPU 可達 4 GHz),接下來窮舉速度就不堪負荷。

時間複雜度 $O(N!)$ , 空間複雜度 $O(N)$

踏出第一步-動態規劃

當學會動態規劃後,解決這一 NP-Complete 問題時,我們使用記憶體空間去掉重複換取時間,定義狀態 dp[i][j] 表示當前走過的點集合 $i$ ,停留最後一個點位置 $j$ ,點集合可以用 bitmask 的技巧壓縮,集合狀態個數 $2^N$ ,窮舉所有停留點轉移需要 $O(N)$ 次。當 $N \ge 25$ 時,理論次數達 8 億之多,再加上記憶體存取效能,導致運行約數十秒。

時間複雜度 $O(N^2 \cdot 2^N)$ , 空間複雜度 $O(N \cdot 2^N)$

更進一步

我們發現到基礎窮舉中可以剪枝卻又不夠局觀,在動態規劃中擁有移除重複卻又不能剪枝,動態規劃又佔據大量的記憶體空間,那麼使用更厲害的 branch-and-reduce 算法吧。簡單介紹一下 branch-and-reduce 算法的特性:

  • 屬於窮舉的一環
  • 在每一步中,花費較多的時間估算整體基礎花費,用來剪枝
  • 窮舉一步之後,將盤面選擇縮小,甚至做到同構移除來減少窮舉分支,這個步驟稱為 reduce
  • 使用的記憶體空間少 (與一般窮舉相當)

運用在 TSP 問題時,我們首先將輸入轉換成 $N \times N$ 矩陣,當移動從 $u$$v$ 時,矩陣就會變成 $(N-1) \times (N-1)$ 大小。

計算基礎花費

對於一個 TSP 而言,每一個點必然有入邊和出邊,因此基礎花費必然每一個節點權重最小的入邊和出邊總和。範例操作如下:

$$\begin{align*} \text{cost matrix } M_5 = \begin{bmatrix} \infty & 20 & 30 & 10 & 11 \\ 15 & \infty & 16 & 4 & 2 \\ 3 & 5 & \infty & 2 & 4 \\ 19 & 6 & 18 & \infty & 3 \\ 16 & 4 & 7 & 16 & \infty \end{bmatrix} \end{align*}$$

計算出邊最小值,將每一行的最小值算出,例如第一行 $[\infty \; 20 \; 30 \; 10 \; 11]$ 為 10,最後分別扣掉 10, 2, 2, 3, 4,得到出邊基礎花費 $10+2+2+3+4 = 21$ ,接著我們扣除掉基礎花費得到

$$\begin{align*} \text{reduce row, cost matrix } M_5 = \begin{bmatrix} \infty & 10 & 20 & 0 & 1 \\ 13 & \infty & 14 & 2 & 0 \\ 1 & 3 & \infty & 0 & 2 \\ 16 & 3 & 15 & \infty & 0 \\ 12 & 0 & 3 & 12 & \infty \end{bmatrix} \end{align*}$$

此時,我們再針對入邊找到基礎花費,將每一個列的最小值算出,分別為 1, 0, 3, 0, 0,因此得到基礎花費 $21+1+3 = 25$ ,如此一來當下限高於當前的最佳解就可以退出。下一階段的矩陣就會被改造成如下所示

$$\begin{align*} \text{reduce row/column, cost matrix }M_5 = \begin{bmatrix} \infty & 10 & 17 & 0 & 1 \\ 12 & \infty & 11 & 2 & 0 \\ 0 & 3 & \infty & 0 & 2 \\ 15 & 3 & 12 & \infty & 0 \\ 11 & 0 & 0 & 12 & \infty \end{bmatrix} \end{align*}$$

當我們從某一點 $u$ 移動到 $v$ 時,移除掉矩陣中的某一行 $u$ 和某一列 $v$ 即可。特別注意到,上述範例為完全圖,對於稀疏圖會發生無解情況,需要特別判斷。

優化一夜

用陣列實作不算困難,當抽出一行一列時,可能會需要複製整個矩陣,又或者只用標記來防止計算到經過的點,不管哪一種寫法都不算漂亮。這裡提供一個解決方法:在處理基礎花費時,對整個矩陣每行每列扣一個定值,實作時不真正修改數值,而是把差值記錄下來,用 $O(N)$ 空間 $Dx[], \; Dy[]$ 記錄差值,如此一來就不用耗費 $O(N^2)$ 空間和時間進行複製。

相對地,這會造成計算真正值 $M'_{x, y} = M_{x, y} - Dx[x] - Dy[y]$ 需要數次的記憶體存取。無妨地,我們已經將了每一層的空間和搬運次數,又因為陣列小到可以全塞進快取中。

優化二夜

對於稀疏圖而言,用矩陣實作過於浪費空間,而且很容易拖累速度,因此可以使用雙十字鏈表 Dancing Links 來模擬上述的矩陣,高效地降低運行常數。實作複雜度指數上升。

優化三夜

遞迴窮舉必然會遇到復原操作,除非我們完整地複製每一層狀態。通常我們只會修改差異,這裡需要對整個矩陣修改,這導致差異數量非常多。為了加速復原操作,遞迴時維護自己 stack 來復原修改差異,而且僅當 $Dx[], Dy[]$ 不為零的時候才進行儲存,因為一部分的 reduce 結果會是零,這些影響的效能非常多。

除了上述外,計算基礎花費仍佔據了 $O(N^2)$ 的時間並佔大部分的運行時間。為了減少計算量,可以做到以下幾點:

  • 在每一行/列計算最小值時,走訪 Dancing Links 每一行每一列,如果當前最小值已經被更新到 0 時,可以直接退出。
  • 傳遞當前已知最佳解,當累積計算花費總和高於已知解就退出。

優化四夜

從啟發的觀點來看,每一階段有 $n$ 種選擇,各別算出預期最少花費後,排序後優先從花費小的開始拓展。

下述程式只支援頂點個數 30 以內,在頂點個數 30 時,只需要耗費 161 ms 即可完成。
當我將頂點個數放大 50 時,小夥伴們的程序花了 2 個小時還沒結果,而我寫的版本跑了 1X 秒就結束。測資的極限就只到這了,暫時先擱著吧。題目鏈結 20007. Fast Travelling Salesman Problem

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
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
#include <bits/stdc++.h>
using namespace std;

#define MAXNODE 5000
#define MAXCOL 32
#define MAXN 32
#define INF 0x3f3f3f
struct DancingNode {
int left, right, up, down;
int ch, rh, w;
} DL[MAXNODE];
struct HelperNode {
int head, size, next, prev;
} HN[MAXNODE];
int help_head;
int s[MAXCOL];
int head, size, Ans;
void Remove(int c) {
DL[DL[c].right].left = DL[c].left;
DL[DL[c].left].right = DL[c].right;
for (int i = DL[c].down; i != c; i = DL[i].down) {
if (HN[DL[i].rh].head == i)
HN[DL[i].rh].head = DL[i].left;
HN[DL[i].rh].size--;
DL[DL[i].right].left = DL[i].left;
DL[DL[i].left].right = DL[i].right;
s[DL[i].ch]--;
}
}
void Resume(int c) {
for (int i = DL[c].down; i != c; i = DL[i].down) {
if (HN[DL[i].rh].head == i)
HN[DL[i].rh].head = DL[i].right;
HN[DL[i].rh].size++;
DL[DL[i].right].left = i;
DL[DL[i].left].right = i;
s[DL[i].ch]++;
}
DL[DL[c].right].left = c;
DL[DL[c].left].right = c;
}
void Reduce(int i) {
int t = DL[i].rh;
HN[HN[t].next].prev = HN[t].prev;
HN[HN[t].prev].next = HN[t].next;
s[DL[i].ch]--;

DL[DL[i].down].up = DL[i].up;
DL[DL[i].up].down = DL[i].down;
for (int k = DL[i].right; k != i; k = DL[k].right) {
DL[DL[k].down].up = DL[k].up;
DL[DL[k].up].down = DL[k].down;
s[DL[k].ch]--;
}
}
void Recover(int i) {
int t = DL[i].rh;
HN[HN[t].next].prev = t;
HN[HN[t].prev].next = t;
s[DL[i].ch]++;

DL[DL[i].down].up = i;
DL[DL[i].up].down = i;
for (int k = DL[i].right; k != i; k = DL[k].right) {
DL[DL[k].down].up = k;
DL[DL[k].up].down = k;
s[DL[k].ch]++;
}
}

void Select(int i, int j) {
int s = DL[i].rh;
Reduce(i);
Remove(j);
HN[HN[s].next].prev = HN[s].prev;
HN[HN[s].prev].next = HN[s].next;
}

void Cancel(int i, int j) {
int s = DL[i].rh;
Resume(j);
Recover(i);
HN[HN[s].next].prev = s;
HN[HN[s].prev].next = s;
}

int new_node(int up, int down, int left, int right) {
DL[size].up = up, DL[size].down = down;
DL[size].left = left, DL[size].right = right;
DL[up].down = DL[down].up = DL[left].right = DL[right].left = size;
return size++;
}

void new_row(int n, int R[][2], int rh) {
int r, row = -1, k;
int h = size;
for (int i = 0; i < n; i++) {
r = R[i][0];
DL[size].ch = r, s[r]++;
DL[size].rh = rh;
DL[size].w = R[i][1];
if (row == -1) {
row = new_node(DL[DL[r].ch].up, DL[r].ch, size, size);
} else {
k = new_node(DL[DL[r].ch].up, DL[r].ch, DL[row].left, row);
}
}
HN[rh].size = n;
HN[rh].head = h;
HN[rh].next = HN[help_head].next;
HN[rh].prev = help_head;
HN[HN[help_head].next].prev = rh;
HN[help_head].next = rh;
}
void init(int m) {
size = 0;
help_head = 0, HN[help_head].next = HN[help_head].prev = help_head;
head = new_node(0, 0, 0, 0);
for (int i = 1; i <= m; i++) {
new_node(i, i, DL[head].left, head);
DL[i].ch = i, s[i] = 0;
DL[i].rh = 0; // important pointer (fail pointer)
}
}

int Dx[MAXN], Dy[MAXN];
int markRStk[MAXNODE], markRidx = -1;
int markCStk[MAXNODE], markCidx = -1;

void printDL(int n) {
for (int i = HN[help_head].prev; i != help_head; i = HN[i].prev) {
int j = HN[i].head, prev = 1;
printf("%2d ", DL[j].rh);
do {
while (prev < DL[j].ch)
prev++, printf(" -- ");
prev = DL[j].ch+1;
int rh = DL[j].rh, ch = DL[j].ch;
printf("[%2d]", DL[j].w - Dx[rh] - Dy[ch]);
j = DL[j].right;
} while (j != HN[i].head);
while (prev <= n)
prev++, printf(" -- ");
puts("");
}
puts("----------------------");
}
int cH(int limit = INF) {
int ins = 0;
for (int i = HN[help_head].prev; i != help_head; i = HN[i].prev) {
if (HN[i].size == 0)
return INF;
int j = HN[i].head, v = INF;
int rh = DL[j].rh, ch;
do {
ch = DL[j].ch;
v = min(v, DL[j].w - Dx[rh] - Dy[ch]);
j = DL[j].right;
} while (j != HN[i].head && v);
if (v == INF || ins+v >= limit)
return INF;
if (v) {
ins += v;
Dx[rh] += v;
markRStk[++markRidx] = v;
markRStk[++markRidx] = rh;
}
}
/*
printf("cH first\n");
printDL(4);
*/

for (int i = DL[head].right; i != head; i = DL[i].right) {
if (DL[i].down == i)
return INF;
int j = DL[i].down, v = INF;
int ch = DL[i].ch, rh;
do {
rh = DL[j].rh;
v = min(v, DL[j].w - Dx[rh] - Dy[ch]);
j = DL[j].down;
} while (j != i && v);
if (v == INF || ins+v >= limit)
return INF;
if (v) {
ins += v;
Dy[ch] += v;
markCStk[++markCidx] = v;
markCStk[++markCidx] = ch;
}
}
return ins;
}
void rH(int backRidx, int backCidx) {
while (markRidx > backRidx) {
int a = markRStk[markRidx--];
int b = markRStk[markRidx--];
Dx[a] -= b;
}
while (markCidx > backCidx) {
int a = markCStk[markCidx--];
int b = markCStk[markCidx--];
Dy[a] -= b;
}
}

#ifdef SAME_CUT
int isSame(int vx, int vy) {
int i = HN[vx].head, j = HN[vy].head;
do {
if (DL[i].ch != DL[j].ch && (DL[i].ch != vy && DL[j].ch != vx))
return 0;
int e1 = DL[i].w - Dx[vx] - Dy[DL[i].ch];
int e2 = DL[j].w - Dx[vy] - Dy[DL[j].ch];
if (e1 != e2)
return 0;
i = DL[i].right;
j = DL[j].right;
} while (i != HN[vx].head && j != HN[vy].head);
return i == HN[vx].head && j == HN[vy].head;
}
#endif

#ifdef MST_CUT
int parent[MAXN];
int findp(int x) {
return x == parent[x] ? x : (parent[x]=findp(parent[x]));
}
int joint(int x, int y) {
x = findp(x), y = findp(y);
if (x == y) return 0;
parent[y] = x;
return 1;
}
int MST(int st, int ed, int limit = INF) {
static int used[MAXN] = {}, cases = 0;
if (HN[HN[help_head].prev].prev == help_head)
return 0;
int i = HN[help_head].prev;
int n = 0;
int ret = 0, in = INF, out = INF;
cases++;
do {
int j = HN[i].head;
do {
int x = DL[j].rh, y = DL[j].ch;
int w = DL[j].w - Dx[x] - Dy[y];
if (x == st) {
in = min(in, w);
} else if (y == ed) {
out = min(out, w);
} else {
if (used[x] != cases)
used[x] = cases, parent[x] = x, n++;
if (used[y] != cases)
used[y] = cases, parent[y] = y, n++;
n -= joint(x, y);
}
j = DL[j].right;
} while (j != HN[i].head);
i = HN[i].prev;
} while (i != help_head);
ret = in + out;
if (ret >= limit)
return limit;
if (n > 1) {
fprintf(stderr, "cc cut\n");
return INF;
}
return ret;
}
#endif
void DFS(int u, int lower_bound) {
#ifdef MST_CUT
if (lower_bound + MST(u, 1, Ans - lower_bound) >= Ans) {
return ;
}
#else
if (lower_bound >= Ans) {
return ;
}
#endif

if (HN[HN[help_head].prev].prev == help_head) {
Ans = lower_bound;
return ;
}

int i = HN[u].head;
vector< pair<int, int> > pick;
do {
int v = DL[i].ch;
int runnable = 1;
if (v == 1)
runnable = 0;
#ifdef SAME_CUT
if (runnable) {
int e1 = DL[i].w - Dx[u] - Dy[v];
for (int j = HN[u].head; j != i && runnable; j = DL[j].right) {
int e2 = DL[j].w - Dx[u] - Dy[DL[j].ch];
if (e1 != e2)
continue;
if (isSame(v, DL[j].ch))
runnable = 0;
}
}
#endif

if (runnable) {
int backRidx = markRidx, backCidx = markCidx;
int tu = HN[u].head;
Select(tu, v);
int c = cH(Ans - lower_bound);
if (c != INF) {
pick.push_back(make_pair(lower_bound + c + DL[i].w-Dx[u]-Dy[v], v));
// DFS(v, lower_bound + c + DL[i].w-Dx[u]-Dy[v]);
}
rH(backRidx, backCidx);
Cancel(tu, v);
}
i = DL[i].right;
} while (i != HN[u].head);

sort(pick.begin(), pick.end());

for (int i = 0; i < pick.size(); i++) {
int v = pick[i].second;
/*
printf("%d -> %d, %d\n", u, v, pick[i].first);
printDL(4);
*/

int backRidx = markRidx, backCidx = markCidx;
int tu = HN[u].head;
Select(tu, v);
/*
printf("Select\n");
printDL(4);
*/

int c = cH(Ans - lower_bound);
/*
printf("Reduce\n");
printDL(4);
*/

if (c != INF)
DFS(v, pick[i].first);
rH(backRidx, backCidx);
Cancel(tu, v);
}
}

int main() {
int G[MAXN][MAXN];
int n;
while (scanf("%d", &n) == 1) {
init(n);
for (int i = 1; i <= n; i++) {
int R[MAXN][2], Rcnt = 0;
for (int j = 1; j <= n; j++) {
assert(scanf("%d", &G[i][j]) == 1);
if (G[i][j])
R[Rcnt][0] = j, R[Rcnt][1] = G[i][j], Rcnt++;
}
new_row(Rcnt, R, i);
}

int backRidx = markRidx, backCidx = markCidx;
Ans = INF;
DFS(1, cH());
rH(backRidx, backCidx);
printf("%d\n", Ans);

}
return 0;
}