淺談旅行推銷員問題 (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

4
0 10 15 20
10 0 35 25
15 35 0 30
20 25 30 0

Sample Output 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

#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;
}